#-------------------------------------------------------------------------------
# Political Depolarization in German Public Opinion, 1980-2010
#
# AUTHORS:     Simon Munzert
#			   Paul Bauer
# INSTITUTION: Department of Politics and Public Administration
#              University of Konstanz, Germany
#			   Institute of Political Science
#			   University of Bern, Switzerland
#
# DATE:        June 2013
#-------------------------------------------------------------------------------

# Remove everything in workspace
rm(list=ls(all=TRUE))

# Load libraries
x <- c("foreign", "reshape", "sna", "arm", "xtable", "lme4", "stringr", "RColorBrewer", "aspace", "mvtnorm", "scatterplot3d")
lapply(x, require, character.only=T)

# Set working directory
setwd("~/Analysis")


################################################################################
### ANALYSIS - PREPARE DATA
################################################################################

load("polarization_data.RData")

years <- c(1980, 1982, 1984, 1986, 1988, 1990, 1991, 1992, 1994, 1996, 1998, 2000, 2002, 2004, 2006, 2008, 2010)
years.names <- paste("y", years, sep="")
df.raw <- rawdata[,c(2,14:(dim(rawdata)[2]))] # pick polarization variables plus year variable
n.years <- length(years)
n.vars <- dim(df.raw)[2]
n.obs <- dim(df.raw)[1]

var.names <- colnames(rawdata[,14:dim(rawdata)[2]])
var.type <- NULL
var.type <- ifelse(str_detect(var.names, "gender."), 1,
				ifelse(str_detect(var.names, "moral."), 2,
					ifelse(str_detect(var.names, "distribution."), 3,
				4)))
var.type # classify items after issue group
table(var.type)
var.order <- order(var.type)


################################################################################
### DESCRIPTIVE STATISTICS
################################################################################

### ISSUE BY YEAR SUMMARY TABLE
by.year.table <- matrix(NA, n.vars, n.years)
for (var in 2:n.vars){
	a <-table(df.raw$year, df.raw[,var])
	for (y in 1:n.years){
		if (sum(a[y,])>0)
		by.year.table[var,y]<- 1
		else by.year.table[var,y]<- 0
	}
}
#write.table(by.year.table, 'by.year.table.txt')



################################################################################
### ANALYSIS - MULTILEVEL MODEL ISSUE*ISSUE
################################################################################


### BUILD FUNCTION gen.correlation.matrix TO GENERATE CORRELATION MATRIX FOR ANALYSIS
################################################################################

gen.correlation.matrix <- function(dataset.raw, type=c("raw","rich","halfrich","mean")) {
###--- GENERATE BASE VARIABLES ---###
years <- sort(unique(dataset.raw$year))
years.names <- paste("y", years, sep="")
df.vars <- dataset.raw
n.years <- length(years)
n.vars <- dim(df.vars)[2]
n.obs <- dim(df.vars)[1]

###--- MARK ISSUE GROUPS ---###
	# gender issues == 2
	# traditional issues == 2
	# socio-economic issues == 3
	# immigration issues == 4
var.names <- colnames(df.vars[,2:dim(df.vars)[2]])
var.type <- NULL
var.type <- ifelse(str_detect(var.names, "gender."), 1,
				ifelse(str_detect(var.names, "moral."), 2,
					ifelse(str_detect(var.names, "distribution."), 3,
				4)))
var.type # classify items after issue group
table(var.type)
var.order <- order(var.type)

###--- GENERATE CORRELATION MATRIX (ISSUE*ISSUE*YEAR) ---###
n.vars.sub <- n.vars - 1 # n.vars - year variable
allbus.cor <- array (NA, dim =c(n.vars.sub, n.vars.sub, n.years))
for (y in 1:n.years){
	oneyr <- df.vars[df.vars$year==years[y],2:dim(df.vars)[2]]
	for (v1 in 1:n.vars.sub){
		a <- v1+1
		if (a <= n.vars.sub){
		for (v2 in a:n.vars.sub){
			allbus.cor[v1,v2,y] <- cor(oneyr[,v1], oneyr[,v2], use = "pairwise.complete.obs", method="pearson") # issue-issue correlation (Pearson) within one year
			}
		}
 }
 allbus.cor[,,y] <- symmetrize(allbus.cor[,,y], rule='upper') # Copy the upper triangle over the lower triangle
}

###--- GENERATE MEAN CORRELATION MATRIX (ISSUE*ISSUE) ---###
allbus.cor.mean<-matrix(NA, n.vars.sub, n.vars.sub)
for (v1 in var.order){
	for (v2 in var.order){
		allbus.cor.mean[v1,v2] <- mean(allbus.cor[v1,v2,], na.rm=T)
		}
	}
allbus.cor.mean <- symmetrize(allbus.cor.mean, rule='upper')

###--- RESHAPE CORRELATION MATRIX ---###
n.pairs<- n.vars.sub*(n.vars.sub-1)/2
rest.names <- c('v1', 'v2', 'pair', 'var.type.v1', 'var.type.v2', 'within.btw.issue', 'eq.diff.types')
allbus.mat <- array (NA, dim =c(n.pairs, n.years+7))
for (y in 1:n.years){
	pa<-1
	for (v1 in 1:n.vars.sub) {
		for (v2 in 1:n.vars.sub) {
			if (v2>v1){
			allbus.mat[pa,n.years+1] <- v1 # variable number v1
			allbus.mat[pa,n.years+2] <- v2 # variable number v2
			allbus.mat[pa,n.years+3] <- pa # pair indicator
			allbus.mat[pa,n.years+4] <- var.type[v1] # issue group v1
			allbus.mat[pa,n.years+5] <- var.type[v2] # issue group v2
			allbus.mat[pa,n.years+6] <-	ifelse (var.type[v1]==var.type[v2],1,0) # issues in the same group?
			if (var.type[v1]==var.type[v2])			allbus.mat[pa,n.years+7] <- var.type[v1] # pair group indicator
			else if (var.type[v1]==1&var.type[v2]==2)	allbus.mat[pa,n.years+7] <- 5
			else if (var.type[v1]==1&var.type[v2]==3)	allbus.mat[pa,n.years+7] <- 6
			else if (var.type[v1]==1&var.type[v2]==4)	allbus.mat[pa,n.years+7] <- 7
			else if (var.type[v1]==2&var.type[v2]==1)	allbus.mat[pa,n.years+7] <- 8
			else if (var.type[v1]==2&var.type[v2]==3)	allbus.mat[pa,n.years+7] <- 9
			else if (var.type[v1]==2&var.type[v2]==4)	allbus.mat[pa,n.years+7] <- 10
			else if (var.type[v1]==3&var.type[v2]==1)	allbus.mat[pa,n.years+7] <- 11
			else if (var.type[v1]==3&var.type[v2]==2)	allbus.mat[pa,n.years+7] <- 12
			else if (var.type[v1]==3&var.type[v2]==4)	allbus.mat[pa,n.years+7] <- 13
			else if (var.type[v1]==4&var.type[v2]==1)	allbus.mat[pa,n.years+7] <- 14
			else if (var.type[v1]==4&var.type[v2]==2)	allbus.mat[pa,n.years+7] <- 15
			else if (var.type[v1]==4&var.type[v2]==3)	allbus.mat[pa,n.years+7] <- 16
			allbus.mat[pa,y]<-allbus.cor[v1,v2,y]
			pa<-pa+1
				}
			}
		}
	}
label.names <- c(years.names,rest.names)
colnames(allbus.mat) <- label.names
allbus.mat	<- as.data.frame(allbus.mat)

###--- MANUALLY RESHAPE DATA FRAME, LONG ---###
allbus.mat.long <- array (NA, dim=c (8000, 9)) # 8000 is an arbitrarily high number - data.frame is reduced later
r<-1
for (pa in 1:n.pairs){
	for (y in 1:n.years){
		if (!is.na(allbus.mat[pa,y])){
			allbus.mat.long[r,9]<- allbus.mat[pa,y]
			allbus.mat.long[r,8]<- allbus.mat[pa,"eq.diff.types"]
			allbus.mat.long[r,7]<- allbus.mat[pa,"within.btw.issue"]
			allbus.mat.long[r,6]<- y
			allbus.mat.long[r,5]<- years[y]
			allbus.mat.long[r,4]<- allbus.mat[pa,"v2"]
			allbus.mat.long[r,3]<- allbus.mat[pa,"v1"]
			allbus.mat.long[r,2]<- allbus.mat[pa,"pair"]
			allbus.mat.long[r,1]<- r
			r<-r+1
			}
		}
	}
	
allbus.mat.long <- allbus.mat.long[complete.cases(allbus.mat.long),]	
label.names <- c('ones', 'pair', 'v1', 'v2', 'year', 'n.year', 'within.btw.issue', 'eq.diff.types', 'coeff')
colnames(allbus.mat.long) <- label.names
allbus.mat.long	<- as.data.frame(allbus.mat.long)
if (type=="rich") {
return(allbus.mat.long) }
else if (type=="halfrich") {
return(allbus.mat)
}
else if (type=="raw") {
return(allbus.cor)
}
else if (type=="mean") {
return(allbus.cor.mean)
}
else
cat("Error")
}
################################################################################

#####PLOT AVERAGE CORRELATION COEFFICIENTS
################################################################################

matrix.cor.mean <- gen.correlation.matrix(df.raw, type="mean")

### BUILD FUNCTION triangleplot2a TO GENERATE TRIANGLEPLOTS
################################################################################

# triangleplot2a is a slight modification of the triangleplot() function of the package 'arm'
triangleplot2a <- function (x, y = NULL, cutpts = NULL, details = TRUE, n.col.legend = 5, 
    cex.col = 0.7, cex.var = 0.9, digits = 1, labels.cutpts= NULL, color = FALSE) 
{
    if (!is.matrix(x)) 
        stop("x must be a matrix!")
    if (dim(x)[1] != dim(x)[2]) 
        stop("x must be a square matrix!")
    x.na <- x
    x.na[is.na(x.na)] <- -999
    z.plot <- x
    if (is.null(y)) {
        z.names <- dimnames(x)[[2]]
    }
    else {
        z.names <- y
    }
    for (i in 1:dim(z.plot)[1]) for (j in i:dim(z.plot)[2]) z.plot[i, 
        j] <- NA
    layout(matrix(c(2, 1), 1, 2, byrow = FALSE), c(10.5, 1.5))
    layout(matrix(c(2, 1), 1, 2, byrow = FALSE), c(10.5, 1.5))
    if (is.null(cutpts)) {
        if (details) {
            neg.check <- abs(sum(z.plot[z.plot < 0], na.rm = T))
            if (neg.check > 0) {
                z.breaks <- sort(c(0, seq(min(z.plot, na.rm = T), 
                  max(z.plot, na.rm = T), length = n.col.legend)))
            }
            else {
                z.breaks <- seq(min(z.plot, na.rm = T), max(z.plot, 
                  na.rm = T), length = n.col.legend + 1)
            }
            for (i in 1:4) {
                n1 <- length(unique(round(z.breaks, digits = digits)))
                n2 <- length(z.breaks)
                ifelse((n1 != n2), digits <- digits + 1, digits <- digits)
            }
            if (digits > 3) {
                stop("Too many digits! Try to adjust n.col.legend to get better presentation!")
            }
        }
        else {
            postive.z <- na.exclude(unique(round(z.plot[z.plot > 
                0], digits = digits)))
            neg.check <- abs(sum(z.plot[z.plot < 0], na.rm = T))
            ifelse(neg.check > 0, negative.z <- na.exclude(unique(round(z.plot[z.plot < 
                0], digits = digits))), negative.z <- 0)
            max.z <- max(z.plot, na.rm = T)
            min.z <- min(z.plot, na.rm = T)
            z.breaks <- sort(unique(c(postive.z, negative.z)))
            n.breaks <- length(z.breaks)
            l.legend <- ceiling(n.col.legend/2)
            if (n.breaks > 8) {
                if (neg.check > 0) {
                  postive.z <- seq(0, max(postive.z), length = l.legend + 
                    1)
                  negative.z <- seq(min(negative.z), 0, length = l.legend)
                  z.breaks <- sort(unique(c(postive.z, negative.z)))
                  n.breaks <- length(z.breaks)
                  z.breaks[1] <- min.z
                  z.breaks[n.breaks] <- max.z
                  n.col.legend <- length(z.breaks) - 1
                }
                else {
                  postive.z <- seq(0, max(postive.z), length = n.col.legend + 
                    1)
                  z.breaks <- sort(unique(c(postive.z, negative.z)))
                  n.breaks <- length(z.breaks)
                  z.breaks[1] <- min.z
                  z.breaks[n.breaks] <- max.z
                  n.col.legend <- length(z.breaks) - 1
                }
            }
            else {
                if (neg.check > 0) {
                  z.breaks <- sort(c(0, seq(min(z.plot, na.rm = T), 
                    max(z.plot, na.rm = T), length = n.col.legend)))
                }
                else {
                  z.breaks <- seq(min(z.plot, na.rm = T), max(z.plot, 
                    na.rm = T), length = n.col.legend + 1)
                }
            }
        }
    }
    if (!is.null(cutpts)) {
        z.breaks = cutpts
        n.breaks <- length(z.breaks)
        n.col.legend <- length(z.breaks) - 1
    }
    if (color) {
        z.colors <- brewer.pal(n.breaks-1, "RdYlGn")
		 
    }
    else {
        z.colors <- gray(n.col.legend:1/n.col.legend)
    }
    par(mar = c(0.5, 0.1, 2, 0.1), pty = "m")
    plot(c(0, 1), c(min(z.breaks), max(z.breaks)), type = "n", 
        bty = "n", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
    for (i in 2:(length(z.breaks))) {
        rect(xleft = 0.5, ybottom = z.breaks[i - 1], xright = 1, 
            ytop = z.breaks[i], col = z.colors[i - 1])
        text(x = 0.45, y = z.breaks[i - 1], labels = format(round(z.breaks[i - 
            1], digits)), cex = cex.col, adj = 1, xpd = TRUE)
    }
    rect(xleft = 0.5, ybottom = z.breaks[length(z.breaks)], xright = 1, 
        ytop = z.breaks[length(z.breaks)], col = z.colors[length(z.colors)])
    text(x = 0.45, y = z.breaks[length(z.breaks)], labels = format(round(z.breaks[length(z.breaks)], 
        digits)), cex = cex.col, adj = 1, xpd = TRUE)
    par(mar = c(0.1, 0.1, 2, 0.1), pty = "m")
    image(x = 1:dim(z.plot)[1], y = 1:dim(z.plot)[2], z = z.plot, 
        xaxt = "n", yaxt = "n", bty = "n", col = z.colors, breaks = z.breaks, 
        xlim = c(-2, dim(z.plot)[1] + 0.5), ylim = c(-1, dim(z.plot)[2] + 
            0.5), xlab = "", ylab = "")
    text(x = 1:dim(z.plot)[1], y = 1:dim(z.plot)[2], labels = z.names, 
        cex = cex.var, adj = 1, xpd = TRUE)
    for (i in 1:dim(z.plot)[1]) {
        for (j in i:dim(z.plot)[2]) {
            if (x.na[i, j] == -999 & i != j) 
                points(x = j, y = i, pch = "x", cex = 0.9)
        }
    }
    for (i in 1:dim(z.plot)[1]) {
        for (j in i:dim(z.plot)[2]) {
            if (x.na[i, j] != -999 & i != j & x.na[i, j] > 0) 
                points(x = j, y = i, pch = "+", cex = 0.9)
        }
    }
    for (i in 1:dim(z.plot)[1]) {
        for (j in i:dim(z.plot)[2]) {
            if (x.na[i, j] != -999 & i != j & x.na[i, j] < 0) 
                points(x = j, y = i, pch = "-", cex = 0.9)
        }
    }
}


################################################################################
### CORRELATION PLOT FOR ALL ATTITUDE SCALE PAIRS (FIGURE 1)
################################################################################

# pdf(file="triangleplot2a_corr.pdf", height=8, width=11, family="URWTimes")
par(oma=c(0,1,0,0))
par (mar=c(1.5,1.5,0,1.5))
par(mfrow=c(1,1))
triangleplot2a(matrix.cor.mean, var.names[var.order], color=T, cex.var=.9, n.col.legend=10,  cutpts=c(-.1,-.1,0,0,0,0,.1,.2,.3,.4,.5,.6))
arrows(1.5,0.5,6.5,0.5,col="black",length=0, lwd=2) # arrows for gender items
arrows(6.5,0.5,6.5,5.5,col="black",length=0, lwd=2)
	for (i in seq(1.5,5.5,1)) {
	arrows(i,i-1,i,i,col="black",length=0, lwd=2)
	arrows(i,i,i+1,i,col="black",length=0, lwd=2)
	}
arrows(7.5,6.5,11.5,6.5,col="black",length=0, lwd=2)
arrows(11.5,6.5,11.5,10.5,col="black",length=0, lwd=2)
	for (i in seq(7.5,10.5,1)) {
	arrows(i,i-1,i,i,col="black",length=0, lwd=2)
	arrows(i,i,i+1,i,col="black",length=0, lwd=2)
	}
arrows(12.5,11.5,18.5,11.5,col="black",length=0, lwd=2)
arrows(18.5,11.5,18.5,17.5,col="black",length=0, lwd=2)
	for (i in seq(12.5,17.5,1)) {
	arrows(i,i-1,i,i,col="black",length=0, lwd=2)
	arrows(i,i,i+1,i,col="black",length=0, lwd=2)
	}
arrows(19.5,18.5,24.5,18.5,col="black",length=0, lwd=2)
arrows(24.5,18.5,24.5,23.5,col="black",length=0, lwd=2)
	for (i in seq(19.5,23.5,1)) {
	arrows(i,i-1,i,i,col="black",length=0, lwd=2)
	arrows(i,i,i+1,i,col="black",length=0, lwd=2)
	}
# dev.off()



################################################################################
#####MULTILEVEL MODELS ISSUE*ISSUE
################################################################################

### prepare dataset for analysis
df.analysis <- gen.correlation.matrix(df.raw, type="rich")

issue.types.within <- df.analysis$within.btw.issue
issue.types.all <- df.analysis$eq.diff.types
leng <- length(df.analysis$ones)
issue.types.equal <- vector(mode='numeric', leng) # 4 different same-issue-group pairs + mixed-issue pairs (5)
for (i in 1:leng){
	issue.types.equal[i] <- ifelse(df.analysis$eq.diff.types[i] <= 4, df.analysis$eq.diff.types[i], 5)
	}
rho <- df.analysis$coeff
t <- (df.analysis$year - 1994) / 10 #centered at 1994 and by 10 years
pair <- df.analysis$pair
df.models <- data.frame(cbind(rho, t, pair, issue.types.within, issue.types.equal)) # create model data frame


### Model results
(M1 <- lmer (rho ~ 1 + t + (1+t|pair)))
(M1.poly2 <- lmer (rho ~ 1 + poly(t, 2, raw=T) + (1+t|pair)))
(M1.poly3 <- lmer (rho ~ 1 + poly(t, 3, raw=T) + (1+t|pair)))

(M2 <- lmer (rho ~ t*issue.types.within + (t|pair)))

(M3 <- lmer (rho ~ t*factor(issue.types.equal) + (t|pair)))
(M3.poly2 <- lmer (rho ~ poly(t, 2, raw=T)*factor(issue.types.equal) + (t|pair)))
(M3.poly3 <- lmer (rho ~ poly(t, 3, raw=T)*factor(issue.types.equal) + (t|pair)))

### Model results, without intercept and main time effect
(M3a <- lmer (rho ~ t*factor(issue.types.equal) + (t|pair) -1 -t))
(M3a.poly2 <- lmer (rho ~ poly(t, 2, raw=T)*factor(issue.types.equal) + (t|pair) -1 -t))
(M3a.poly3 <- lmer (rho ~ poly(t, 3, raw=T)*factor(issue.types.equal) + (t|pair) -1 -t))




################################################################################
### EFFECT PLOT FOR ISSUE DOMAIN POLARIZATION TRENDS (FIGURE 3)
################################################################################

M3a.women.hi <- vector()
M3a.women.lo <- vector()
M3a.moral.hi <- vector()
M3a.moral.lo <- vector()
M3a.distribution.hi <- vector()
M3a.distribution.lo <- vector()
M3a.immigration.hi <- vector()
M3a.immigration.lo <- vector()
for (i in 1:length(years)) {
M3a.women.hi[i] <- (M3a@fixef[1])+(M3a@fixef[6])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a)[1])^2+(se.fixef(M3a)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a)[1,6]*(years[i]-1994)/10))
M3a.women.lo[i] <- (M3a@fixef[1])+(M3a@fixef[6])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a)[1])^2+(se.fixef(M3a)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a)[1,6]*(years[i]-1994)/10))
M3a.moral.hi[i] <- (M3a@fixef[2])+(M3a@fixef[7])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a)[2])^2+(se.fixef(M3a)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a)[2,7]*(years[i]-1994)/10))
M3a.moral.lo[i] <- (M3a@fixef[2])+(M3a@fixef[7])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a)[2])^2+(se.fixef(M3a)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a)[2,7]*(years[i]-1994)/10))
M3a.distribution.hi[i] <- (M3a@fixef[3])+(M3a@fixef[8])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a)[3])^2+(se.fixef(M3a)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a)[3,8]*(years[i]-1994)/10))
M3a.distribution.lo[i] <- (M3a@fixef[3])+(M3a@fixef[8])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a)[3])^2+(se.fixef(M3a)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a)[3,8]*(years[i]-1994)/10))
M3a.immigration.hi[i] <- (M3a@fixef[4])+(M3a@fixef[9])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a)[4])^2+(se.fixef(M3a)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a)[4,9]*(years[i]-1994)/10))
M3a.immigration.lo[i] <- (M3a@fixef[4])+(M3a@fixef[9])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a)[4])^2+(se.fixef(M3a)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a)[4,9]*(years[i]-1994)/10))
}

# pdf(file="issue_trends_ci_grey.pdf", height=7, width=7, family="URWTimes")
lightgrey.s <- rgb(180,180,180, 100, max=255)
darkgrey.s <- rgb(50,50,50, 100, max=255)
par(oma=c(1,1,1,0))
par (mar=c(4.5,4.5,2.5,2.5))
par(mfrow=c(1,1))
curve (M3a@fixef[1]+ M3a@fixef[6]*(x-1994)/10,  ylim=c(0.1, .5), xlim=c(1980,2010), lwd=2,  lty='solid', col='white', ylab='Average correlation coefficient', xlab='Year', main='',xaxs = "i", yaxs = "i")
polygon(c(years,rev(years)),c(M3a.women.hi,rev(M3a.women.lo)),col=darkgrey.s, border=NA)
polygon(c(years,rev(years)),c(M3a.moral.hi,rev(M3a.moral.lo)),col=lightgrey.s, border=NA) 
polygon(c(years,rev(years)),c(M3a.distribution.hi,rev(M3a.distribution.lo)),col=darkgrey.s, border=NA) 
polygon(c(years,rev(years)),c(M3a.immigration.hi,rev(M3a.immigration.lo)),col=lightgrey.s, border=NA) 
curve (M3a@fixef[1]+ M3a@fixef[6]*(x-1994)/10,  ylim=c(0.1, .5), xlim=c(1980,2010), lwd=2,  lty='solid', col='white', ylab='Average correlation coefficient', xlab='Year', main='',xaxs = "i", yaxs = "i", add=T)
curve (M3a@fixef[2]+ M3a@fixef[7]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='solid', add=T, xaxs = "i", yaxs = "i")
curve (M3a@fixef[3]+ M3a@fixef[8]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='white', lty='solid', add=T, xaxs = "i", yaxs = "i")
curve (M3a@fixef[4]+ M3a@fixef[9]*(x-1994)/10, xlim=c(1980,2010), lwd=2,  col='black', lty='solid', add=T, xaxs = "i", yaxs = "i")
text(2006,.4,"Gender", srt=16, col="white")
text(1992,.295,"Distribution", srt=-14, col="white")
text(1984,.41,"Immigration", srt=-15, col="black")
text(1998,.217,"Moral", srt=-20, col="black")
# dev.off()


#######################################################################################################
### EFFECT PLOT FOR ISSUE DOMAIN POLARIZATION TRENDS, POLYNOMIAL SPECIFICATION, 3RD DEGREE (FIGURE 7)
#######################################################################################################

# pdf(file="issue_trends_poly3_sims_ci_grey.pdf", height=7, width=7, family="URWTimes")
par(oma=c(1,1,1,0))
par (mar=c(4.5,4.5,2.5,2.5))
par(mfrow=c(1,1))
curve (M3.poly3@fixef[1] + M3.poly3@fixef[2]*((x-1994)/10) + M3.poly3@fixef[3]*((x-1994)/10)^2 + M3.poly3@fixef[4]*((x-1994)/10)^3,  ylim=c(0.1, .5), xlim=c(1980,2010), lwd=2,  lty='solid', col='white', ylab='Average correlation coefficient', xlab='Year', main='',xaxs = "i", yaxs = "i")
polygon(c(years,rev(years)),c(M3.poly3.sim.pred.women.hi,rev(M3.poly3.sim.pred.women.lo)),col=darkgrey.s, border=NA)
polygon(c(years,rev(years)),c(M3.poly3.sim.pred.econ.hi,rev(M3.poly3.sim.pred.econ.lo)),col=darkgrey.s, border=NA) 
polygon(c(years,rev(years)),c(M3.poly3.sim.pred.immig.hi,rev(M3.poly3.sim.pred.immig.lo)),col=lightgrey.s, border=NA) 
curve (M3.poly3@fixef[1] + M3.poly3@fixef[2]*((x-1994)/10) + M3.poly3@fixef[3]*((x-1994)/10)^2 + M3.poly3@fixef[4]*((x-1994)/10)^3,  ylim=c(0.1, .5), xlim=c(1980,2010), lwd=2,  lty='solid', col='white', add=T)
curve ((M3.poly3@fixef[1]+M3.poly3@fixef[6])+(M3.poly3@fixef[2]+M3.poly3@fixef[12])*((x-1994)/10) + (M3.poly3@fixef[3]+M3.poly3@fixef[13])*((x-1994)/10)^2 + (M3.poly3@fixef[4]+M3.poly3@fixef[14])*((x-1994)/10)^3, xlim=c(1980,2010), lwd=2, col='white', lty='solid', add=T,xaxs = "i", yaxs = "i")
curve ((M3.poly3@fixef[1]+M3.poly3@fixef[7])+(M3.poly3@fixef[2]+M3.poly3@fixef[15])*((x-1994)/10) + (M3.poly3@fixef[3]+M3.poly3@fixef[16])*((x-1994)/10)^2 + (M3.poly3@fixef[4]+M3.poly3@fixef[17])*((x-1994)/10)^3, xlim=c(1980,2010), lwd=2,  col='black', lty='solid', add=T,xaxs = "i", yaxs = "i")
text(1994,.35,"Gender", srt=atan_d(M3@fixef[2]*10-.12), col="white")
text(1988,.40,"Immigration", srt=-14, col="black")
text(1995,.29,"Distribution", srt=-18, col="white")
# dev.off()




################################################################################
### CORRELATION TREND PLOT FOR ALL ATTITUDE SCALE PAIRS (FIGURE 2)
################################################################################


df.corr.trend.M1 <- gen.correlation.matrix(df.raw, type="halfrich")
a <- as.matrix(coef(M1)[[1]])
t.pairs <- as.numeric(row.names(a))
num.t.pairs <- length(t.pairs)
for (i in 1:num.t.pairs){
	df.corr.trend.M1[t.pairs[i],n.years+8] <- a[i,2]
	}
df.corr.trend.by.var.M1 <- matrix(NA, n.vars-1, n.vars-1)
for (v1 in 1:n.vars-1){
	for (v2 in 1:n.vars-1){
		if (v1<v2)
		df.corr.trend.by.var.M1[v1,v2] <- df.corr.trend.M1[df.corr.trend.M1[,n.years+1]==v1 & df.corr.trend.M1[,n.years+2]==v2,n.years+8]
		}
	}
df.corr.trend.by.var.M1 <- symmetrize(df.corr.trend.by.var.M1, rule='upper')
df.corr.trend.by.var.M1[df.corr.trend.by.var.M1 <= -.12] <- -.119 # adapt single observation with corr < -.12


# pdf(file="triangleplot2a_trendcorr.pdf", height=8, width=11, family="URWTimes")
par(oma=c(0,1,0,0))
par (mar=c(1.5,1.5,0,1.5))
par(mfrow=c(1,1))
triangleplot2a(df.corr.trend.by.var.M1, var.names[var.order], color=T, cex.var=.9, n.col.legend=10,  cutpts=c(-.12,  -.10, -.08,-.06,-.04,-.02,0, 0,.02,.02, .04, .06), digits=3, labels.cutpts= NULL)
arrows(1.5,0.5,6.5,0.5,col="black",length=0, lwd=2) # arrows for gender items
arrows(6.5,0.5,6.5,5.5,col="black",length=0, lwd=2)
	for (i in seq(1.5,5.5,1)) {
	arrows(i,i-1,i,i,col="black",length=0, lwd=2)
	arrows(i,i,i+1,i,col="black",length=0, lwd=2)
	}
arrows(7.5,6.5,11.5,6.5,col="black",length=0, lwd=2)
arrows(11.5,6.5,11.5,10.5,col="black",length=0, lwd=2)
	for (i in seq(7.5,10.5,1)) {
	arrows(i,i-1,i,i,col="black",length=0, lwd=2)
	arrows(i,i,i+1,i,col="black",length=0, lwd=2)
	}
arrows(12.5,11.5,18.5,11.5,col="black",length=0, lwd=2)
arrows(18.5,11.5,18.5,17.5,col="black",length=0, lwd=2)
	for (i in seq(12.5,17.5,1)) {
	arrows(i,i-1,i,i,col="black",length=0, lwd=2)
	arrows(i,i,i+1,i,col="black",length=0, lwd=2)
	}
arrows(19.5,18.5,24.5,18.5,col="black",length=0, lwd=2)
arrows(24.5,18.5,24.5,23.5,col="black",length=0, lwd=2)
	for (i in seq(19.5,23.5,1)) {
	arrows(i,i-1,i,i,col="black",length=0, lwd=2)
	arrows(i,i,i+1,i,col="black",length=0, lwd=2)
	}
	mtext("<", side=4, at=0.22, line=1.5, las=2, outer=F, cex=.8)
# dev.off()



################################################################################
### ANALYSIS - MULTILEVEL MODEL ISSUE*ISSUE, SUBGROUPS
################################################################################

### generate data subsets: subgroups
df.male <- subset(rawdata, male==1, select=c(2,14:37))
df.female <- subset(rawdata, male==0, select=c(2,14:37))
df.edu.high <- subset(rawdata, education==4 | education==5, select=c(2,14:37))
df.edu.low <- subset(rawdata, education==1 | education==2, select=c(2,14:37))
df.inc.high <- subset(rawdata, income > 2000, select=c(2,14:37))
df.inc.low <- subset(rawdata, income < 500, select=c(2,14:37))
df.int.high <- subset(rawdata, polinterest1 ==1 | polinterest1 == 2, select=c(2,14:37))
df.int.low <- subset(rawdata, polinterest1 == 4 | polinterest1 == 5, select=c(2,14:37))
df.church.pro <- subset(rawdata, confession==1 | confession==2, select=c(2,14:37))
df.church.cat <- subset(rawdata, confession==3, select=c(2,14:37))
df.church.noc <- subset(rawdata, confession==6, select=c(2,14:37))
df.east <- subset(rawdata, west == 0 & year>1990, select=c(2,14:37))
df.west <- subset(rawdata, west == 1, select=c(2,14:37))
df.left <- subset(rawdata, leftright <5, select=c(2,14:37))
df.right <- subset(rawdata, leftright >6, select=c(2,14:37))

### write function for subgroup analyses
subgroup.model <- function(dataframe, model="M1") {
df.analysis <- gen.correlation.matrix(dataframe, type="rich")
issue.types.within <- df.analysis$within.btw.issue
issue.types.all <- df.analysis$eq.diff.types
leng <- length(df.analysis$ones)
issue.types.equal <- vector(mode='numeric', leng) # 4 different same-issue-group pairs + mixed-issue pairs (5)
for (i in 1:leng){
	issue.types.equal[i] <- ifelse(df.analysis$eq.diff.types[i] <= 4, df.analysis$eq.diff.types[i], 5)
	}
rho <- df.analysis$coeff
t <- (df.analysis$year - 1994) / 10 #centered at 1994 and by 10 years
pair <- df.analysis$pair
df.models <- data.frame(cbind(rho, t, pair, issue.types.within, issue.types.equal)) # create model data frame
### Model results
if (model=="M1") {
return(lmer (rho ~ 1 + t + (1+t|pair))) }
else if (model=="M3") {
return(lmer (rho ~ t*factor(issue.types.equal) + (t|pair))) }
else {
return(lmer (rho ~ t*factor(issue.types.equal) + (t|pair) -1 -t))
}
}

### apply function for subgroup analyses
M1.male <- subgroup.model(df.male, model="M1")
M3.male <- subgroup.model(df.male, model="M3")
M3a.male <- subgroup.model(df.male, model="M3a")

M1.female <- subgroup.model(df.female, model="M1")
M3.female <- subgroup.model(df.female, model="M3")
M3a.female <- subgroup.model(df.female, model="M3a")

M1.edu.high <- subgroup.model(df.edu.high, model="M1")
M3.edu.high <- subgroup.model(df.edu.high, model="M3")
M3a.edu.high <- subgroup.model(df.edu.high, model="M3a")

M1.edu.low <- subgroup.model(df.edu.low, model="M1")
M3.edu.low <- subgroup.model(df.edu.low, model="M3")
M3a.edu.low <- subgroup.model(df.edu.low, model="M3a")

M1.inc.high <- subgroup.model(df.inc.high, model="M1")
M3.inc.high <- subgroup.model(df.inc.high, model="M3")
M3a.inc.high <- subgroup.model(df.inc.high, model="M3a")

M1.inc.low <- subgroup.model(df.inc.low, model="M1")
M3.inc.low <- subgroup.model(df.inc.low, model="M3")
M3a.inc.low <- subgroup.model(df.inc.low, model="M3a")

M1.int.high <- subgroup.model(df.int.high, model="M1")
M3.int.high <- subgroup.model(df.int.high, model="M3")
M3a.int.high <- subgroup.model(df.int.high, model="M3a")

M1.int.low <- subgroup.model(df.int.low, model="M1")
M3.int.low <- subgroup.model(df.int.low, model="M3")
M3a.int.low <- subgroup.model(df.int.low, model="M3a")

M1.church.pro <- subgroup.model(df.church.pro, model="M1")
M3.church.pro <- subgroup.model(df.church.pro, model="M3")
M3a.church.pro <- subgroup.model(df.church.pro, model="M3a")

M1.church.cat <- subgroup.model(df.church.cat, model="M1")
M3.church.cat <- subgroup.model(df.church.cat, model="M3")
M3a.church.cat <- subgroup.model(df.church.cat, model="M3a")

M1.church.noc <- subgroup.model(df.church.noc, model="M1")
M3.church.noc <- subgroup.model(df.church.noc, model="M3")
M3a.church.noc <- subgroup.model(df.church.noc, model="M3a")

M1.west <- subgroup.model(df.west, model="M1")
M3.west <- subgroup.model(df.west, model="M3")
M3a.west <- subgroup.model(df.west, model="M3a")

M1.east <- subgroup.model(df.east, model="M1")
M3.east <- subgroup.model(df.east, model="M3")
M3a.east <- subgroup.model(df.east, model="M3a")




################################################################################
### TREND PLOTS FOR SUBPOPULATIONS (FIGURE 4)
################################################################################

white.s <- rgb(230,230,230, 180, max=255)
lightgrey.s <- rgb(209,209,209, 180, max=255)
darkgrey.s <- rgb(140,140,140, 180, max=255)
M1.male.hi <- vector()
M1.male.lo <- vector()
M1.female.hi <- vector()
M1.female.lo <- vector()
M1.edu.high.hi <- vector()
M1.edu.high.lo <- vector()
M1.edu.low.hi <- vector()
M1.edu.low.lo <- vector()
M1.int.high.hi <- vector()
M1.int.high.lo <- vector()
M1.int.low.hi <- vector()
M1.int.low.lo <- vector()
M1.inc.high.hi <- vector()
M1.inc.high.lo <- vector()
M1.inc.low.hi <- vector()
M1.inc.low.lo <- vector()
M1.church.pro.hi <- vector()
M1.church.pro.lo <- vector()
M1.church.cat.hi <- vector()
M1.church.cat.lo <- vector()
M1.church.noc.hi <- vector()
M1.church.noc.lo <- vector()
M1.west.hi <- vector()
M1.west.lo <- vector()
M1.east.hi <- vector()
M1.east.lo <- vector()
for (i in 1:length(years)) {
M1.male.hi[i] <- (M1.male@fixef[1])+(M1.male@fixef[2])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M1.male)[1])^2+(se.fixef(M1.male)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.male)[1,2]*(years[i]-1994)/10))
M1.male.lo[i] <- (M1.male@fixef[1])+(M1.male@fixef[2])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M1.male)[1])^2+(se.fixef(M1.male)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.male)[1,2]*(years[i]-1994)/10))
M1.female.hi[i] <- (M1.female@fixef[1])+(M1.female@fixef[2])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M1.female)[1])^2+(se.fixef(M1.female)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.female)[1,2]*(years[i]-1994)/10))
M1.female.lo[i] <- (M1.female@fixef[1])+(M1.female@fixef[2])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M1.female)[1])^2+(se.fixef(M1.female)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.female)[1,2]*(years[i]-1994)/10))
M1.edu.high.hi[i] <- (M1.edu.high@fixef[1])+(M1.edu.high@fixef[2])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M1.edu.high)[1])^2+(se.fixef(M1.edu.high)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.edu.high)[1,2]*(years[i]-1994)/10))
M1.edu.high.lo[i] <- (M1.edu.high@fixef[1])+(M1.edu.high@fixef[2])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M1.edu.high)[1])^2+(se.fixef(M1.edu.high)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.edu.high)[1,2]*(years[i]-1994)/10))
M1.edu.low.hi[i] <- (M1.edu.low@fixef[1])+(M1.edu.low@fixef[2])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M1.edu.low)[1])^2+(se.fixef(M1.edu.low)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.edu.low)[1,2]*(years[i]-1994)/10))
M1.edu.low.lo[i] <- (M1.edu.low@fixef[1])+(M1.edu.low@fixef[2])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M1.edu.low)[1])^2+(se.fixef(M1.edu.low)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.edu.low)[1,2]*(years[i]-1994)/10))
M1.int.high.hi[i] <- (M1.int.high@fixef[1])+(M1.int.high@fixef[2])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M1.int.high)[1])^2+(se.fixef(M1.int.high)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.int.high)[1,2]*(years[i]-1994)/10))
M1.int.high.lo[i] <- (M1.int.high@fixef[1])+(M1.int.high@fixef[2])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M1.int.high)[1])^2+(se.fixef(M1.int.high)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.int.high)[1,2]*(years[i]-1994)/10))
M1.int.low.hi[i] <- (M1.int.low@fixef[1])+(M1.int.low@fixef[2])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M1.int.low)[1])^2+(se.fixef(M1.int.low)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.int.low)[1,2]*(years[i]-1994)/10))
M1.int.low.lo[i] <- (M1.int.low@fixef[1])+(M1.int.low@fixef[2])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M1.int.low)[1])^2+(se.fixef(M1.int.low)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.int.low)[1,2]*(years[i]-1994)/10))
M1.inc.high.hi[i] <- (M1.inc.high@fixef[1])+(M1.inc.high@fixef[2])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M1.inc.high)[1])^2+(se.fixef(M1.inc.high)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.inc.high)[1,2]*(years[i]-1994)/10))
M1.inc.high.lo[i] <- (M1.inc.high@fixef[1])+(M1.inc.high@fixef[2])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M1.inc.high)[1])^2+(se.fixef(M1.inc.high)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.inc.high)[1,2]*(years[i]-1994)/10))
M1.inc.low.hi[i] <- (M1.inc.low@fixef[1])+(M1.inc.low@fixef[2])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M1.inc.low)[1])^2+(se.fixef(M1.inc.low)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.inc.low)[1,2]*(years[i]-1994)/10))
M1.inc.low.lo[i] <- (M1.inc.low@fixef[1])+(M1.inc.low@fixef[2])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M1.inc.low)[1])^2+(se.fixef(M1.inc.low)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.inc.low)[1,2]*(years[i]-1994)/10))
M1.west.hi[i] <- (M1.west@fixef[1])+(M1.west@fixef[2])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M1.west)[1])^2+(se.fixef(M1.west)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.west)[1,2]*(years[i]-1994)/10))
M1.west.lo[i] <- (M1.west@fixef[1])+(M1.west@fixef[2])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M1.west)[1])^2+(se.fixef(M1.west)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.west)[1,2]*(years[i]-1994)/10))
M1.east.hi[i] <- (M1.east@fixef[1])+(M1.east@fixef[2])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M1.east)[1])^2+(se.fixef(M1.east)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.east)[1,2]*(years[i]-1994)/10))
M1.east.lo[i] <- (M1.east@fixef[1])+(M1.east@fixef[2])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M1.east)[1])^2+(se.fixef(M1.east)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.east)[1,2]*(years[i]-1994)/10))
M1.church.pro.hi[i] <- (M1.church.pro@fixef[1])+(M1.church.pro@fixef[2])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M1.church.pro)[1])^2+(se.fixef(M1.church.pro)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.church.pro)[1,2]*(years[i]-1994)/10))
M1.church.pro.lo[i] <- (M1.church.pro@fixef[1])+(M1.church.pro@fixef[2])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M1.church.pro)[1])^2+(se.fixef(M1.church.pro)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.church.pro)[1,2]*(years[i]-1994)/10))
M1.church.cat.hi[i] <- (M1.church.cat@fixef[1])+(M1.church.cat@fixef[2])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M1.church.cat)[1])^2+(se.fixef(M1.church.cat)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.church.cat)[1,2]*(years[i]-1994)/10))
M1.church.cat.lo[i] <- (M1.church.cat@fixef[1])+(M1.church.cat@fixef[2])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M1.church.cat)[1])^2+(se.fixef(M1.church.cat)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.church.cat)[1,2]*(years[i]-1994)/10))
M1.church.noc.hi[i] <- (M1.church.noc@fixef[1])+(M1.church.noc@fixef[2])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M1.church.noc)[1])^2+(se.fixef(M1.church.noc)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.church.noc)[1,2]*(years[i]-1994)/10))
M1.church.noc.lo[i] <- (M1.church.noc@fixef[1])+(M1.church.noc@fixef[2])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M1.church.noc)[1])^2+(se.fixef(M1.church.noc)[2])^2*((years[i]-1994)/10)^2+2*vcov(M1.church.noc)[1,2]*(years[i]-1994)/10))
}

M3a.female.women.hi <- vector()
M3a.female.women.lo <- vector()
M3a.female.moral.hi <- vector()
M3a.female.moral.lo <- vector()
M3a.female.distribution.hi <- vector()
M3a.female.distribution.lo <- vector()
M3a.female.immigration.hi <- vector()
M3a.female.immigration.lo <- vector()
for (i in 1:length(years)) {
M3a.female.women.hi[i] <- (M3a.female@fixef[1])+(M3a.female@fixef[6])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.female)[1])^2+(se.fixef(M3a.female)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.female)[1,6]*(years[i]-1994)/10))
M3a.female.women.lo[i] <- (M3a.female@fixef[1])+(M3a.female@fixef[6])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.female)[1])^2+(se.fixef(M3a.female)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.female)[1,6]*(years[i]-1994)/10))
M3a.female.moral.hi[i] <- (M3a.female@fixef[2])+(M3a.female@fixef[7])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.female)[2])^2+(se.fixef(M3a.female)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.female)[2,7]*(years[i]-1994)/10))
M3a.female.moral.lo[i] <- (M3a.female@fixef[2])+(M3a.female@fixef[7])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.female)[2])^2+(se.fixef(M3a.female)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.female)[2,7]*(years[i]-1994)/10))
M3a.female.distribution.hi[i] <- (M3a.female@fixef[3])+(M3a.female@fixef[8])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.female)[3])^2+(se.fixef(M3a.female)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.female)[3,8]*(years[i]-1994)/10))
M3a.female.distribution.lo[i] <- (M3a.female@fixef[3])+(M3a.female@fixef[8])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.female)[3])^2+(se.fixef(M3a.female)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.female)[3,8]*(years[i]-1994)/10))
M3a.female.immigration.hi[i] <- (M3a.female@fixef[4])+(M3a.female@fixef[9])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.female)[4])^2+(se.fixef(M3a.female)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.female)[4,9]*(years[i]-1994)/10))
M3a.female.immigration.lo[i] <- (M3a.female@fixef[4])+(M3a.female@fixef[9])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.female)[4])^2+(se.fixef(M3a.female)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.female)[4,9]*(years[i]-1994)/10))
}

M3a.male.women.hi <- vector()
M3a.male.women.lo <- vector()
M3a.male.moral.hi <- vector()
M3a.male.moral.lo <- vector()
M3a.male.distribution.hi <- vector()
M3a.male.distribution.lo <- vector()
M3a.male.immigration.hi <- vector()
M3a.male.immigration.lo <- vector()
for (i in 1:length(years)) {
M3a.male.women.hi[i] <- (M3a.male@fixef[1])+(M3a.male@fixef[6])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.male)[1])^2+(se.fixef(M3a.male)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.male)[1,6]*(years[i]-1994)/10))
M3a.male.women.lo[i] <- (M3a.male@fixef[1])+(M3a.male@fixef[6])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.male)[1])^2+(se.fixef(M3a.male)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.male)[1,6]*(years[i]-1994)/10))
M3a.male.moral.hi[i] <- (M3a.male@fixef[2])+(M3a.male@fixef[7])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.male)[2])^2+(se.fixef(M3a.male)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.male)[2,7]*(years[i]-1994)/10))
M3a.male.moral.lo[i] <- (M3a.male@fixef[2])+(M3a.male@fixef[7])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.male)[2])^2+(se.fixef(M3a.male)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.male)[2,7]*(years[i]-1994)/10))
M3a.male.distribution.hi[i] <- (M3a.male@fixef[3])+(M3a.male@fixef[8])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.male)[3])^2+(se.fixef(M3a.male)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.male)[3,8]*(years[i]-1994)/10))
M3a.male.distribution.lo[i] <- (M3a.male@fixef[3])+(M3a.male@fixef[8])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.male)[3])^2+(se.fixef(M3a.male)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.male)[3,8]*(years[i]-1994)/10))
M3a.male.immigration.hi[i] <- (M3a.male@fixef[4])+(M3a.male@fixef[9])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.male)[4])^2+(se.fixef(M3a.male)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.male)[4,9]*(years[i]-1994)/10))
M3a.male.immigration.lo[i] <- (M3a.male@fixef[4])+(M3a.male@fixef[9])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.male)[4])^2+(se.fixef(M3a.male)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.male)[4,9]*(years[i]-1994)/10))
}

M3a.edu.high.women.hi <- vector()
M3a.edu.high.women.lo <- vector()
M3a.edu.high.moral.hi <- vector()
M3a.edu.high.moral.lo <- vector()
M3a.edu.high.distribution.hi <- vector()
M3a.edu.high.distribution.lo <- vector()
M3a.edu.high.immigration.hi <- vector()
M3a.edu.high.immigration.lo <- vector()
for (i in 1:length(years)) {
M3a.edu.high.women.hi[i] <- (M3a.edu.high@fixef[1])+(M3a.edu.high@fixef[6])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.edu.high)[1])^2+(se.fixef(M3a.edu.high)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.high)[1,6]*(years[i]-1994)/10))
M3a.edu.high.women.lo[i] <- (M3a.edu.high@fixef[1])+(M3a.edu.high@fixef[6])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.edu.high)[1])^2+(se.fixef(M3a.edu.high)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.high)[1,6]*(years[i]-1994)/10))
M3a.edu.high.moral.hi[i] <- (M3a.edu.high@fixef[2])+(M3a.edu.high@fixef[7])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.edu.high)[2])^2+(se.fixef(M3a.edu.high)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.high)[2,7]*(years[i]-1994)/10))
M3a.edu.high.moral.lo[i] <- (M3a.edu.high@fixef[2])+(M3a.edu.high@fixef[7])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.edu.high)[2])^2+(se.fixef(M3a.edu.high)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.high)[2,7]*(years[i]-1994)/10))
M3a.edu.high.distribution.hi[i] <- (M3a.edu.high@fixef[3])+(M3a.edu.high@fixef[8])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.edu.high)[3])^2+(se.fixef(M3a.edu.high)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.high)[3,8]*(years[i]-1994)/10))
M3a.edu.high.distribution.lo[i] <- (M3a.edu.high@fixef[3])+(M3a.edu.high@fixef[8])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.edu.high)[3])^2+(se.fixef(M3a.edu.high)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.high)[3,8]*(years[i]-1994)/10))
M3a.edu.high.immigration.hi[i] <- (M3a.edu.high@fixef[4])+(M3a.edu.high@fixef[9])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.edu.high)[4])^2+(se.fixef(M3a.edu.high)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.high)[4,9]*(years[i]-1994)/10))
M3a.edu.high.immigration.lo[i] <- (M3a.edu.high@fixef[4])+(M3a.edu.high@fixef[9])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.edu.high)[4])^2+(se.fixef(M3a.edu.high)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.high)[4,9]*(years[i]-1994)/10))
}

M3a.edu.low.women.hi <- vector()
M3a.edu.low.women.lo <- vector()
M3a.edu.low.moral.hi <- vector()
M3a.edu.low.moral.lo <- vector()
M3a.edu.low.distribution.hi <- vector()
M3a.edu.low.distribution.lo <- vector()
M3a.edu.low.immigration.hi <- vector()
M3a.edu.low.immigration.lo <- vector()
for (i in 1:length(years)) {
M3a.edu.low.women.hi[i] <- (M3a.edu.low@fixef[1])+(M3a.edu.low@fixef[6])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.edu.low)[1])^2+(se.fixef(M3a.edu.low)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.low)[1,6]*(years[i]-1994)/10))
M3a.edu.low.women.lo[i] <- (M3a.edu.low@fixef[1])+(M3a.edu.low@fixef[6])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.edu.low)[1])^2+(se.fixef(M3a.edu.low)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.low)[1,6]*(years[i]-1994)/10))
M3a.edu.low.moral.hi[i] <- (M3a.edu.low@fixef[2])+(M3a.edu.low@fixef[7])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.edu.low)[2])^2+(se.fixef(M3a.edu.low)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.low)[2,7]*(years[i]-1994)/10))
M3a.edu.low.moral.lo[i] <- (M3a.edu.low@fixef[2])+(M3a.edu.low@fixef[7])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.edu.low)[2])^2+(se.fixef(M3a.edu.low)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.low)[2,7]*(years[i]-1994)/10))
M3a.edu.low.distribution.hi[i] <- (M3a.edu.low@fixef[3])+(M3a.edu.low@fixef[8])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.edu.low)[3])^2+(se.fixef(M3a.edu.low)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.low)[3,8]*(years[i]-1994)/10))
M3a.edu.low.distribution.lo[i] <- (M3a.edu.low@fixef[3])+(M3a.edu.low@fixef[8])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.edu.low)[3])^2+(se.fixef(M3a.edu.low)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.low)[3,8]*(years[i]-1994)/10))
M3a.edu.low.immigration.hi[i] <- (M3a.edu.low@fixef[4])+(M3a.edu.low@fixef[9])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.edu.low)[4])^2+(se.fixef(M3a.edu.low)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.low)[4,9]*(years[i]-1994)/10))
M3a.edu.low.immigration.lo[i] <- (M3a.edu.low@fixef[4])+(M3a.edu.low@fixef[9])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.edu.low)[4])^2+(se.fixef(M3a.edu.low)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.edu.low)[4,9]*(years[i]-1994)/10))
}

M3a.int.high.women.hi <- vector()
M3a.int.high.women.lo <- vector()
M3a.int.high.moral.hi <- vector()
M3a.int.high.moral.lo <- vector()
M3a.int.high.distribution.hi <- vector()
M3a.int.high.distribution.lo <- vector()
M3a.int.high.immigration.hi <- vector()
M3a.int.high.immigration.lo <- vector()
for (i in 1:length(years)) {
M3a.int.high.women.hi[i] <- (M3a.int.high@fixef[1])+(M3a.int.high@fixef[6])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.int.high)[1])^2+(se.fixef(M3a.int.high)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.high)[1,6]*(years[i]-1994)/10))
M3a.int.high.women.lo[i] <- (M3a.int.high@fixef[1])+(M3a.int.high@fixef[6])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.int.high)[1])^2+(se.fixef(M3a.int.high)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.high)[1,6]*(years[i]-1994)/10))
M3a.int.high.moral.hi[i] <- (M3a.int.high@fixef[2])+(M3a.int.high@fixef[7])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.int.high)[2])^2+(se.fixef(M3a.int.high)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.high)[2,7]*(years[i]-1994)/10))
M3a.int.high.moral.lo[i] <- (M3a.int.high@fixef[2])+(M3a.int.high@fixef[7])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.int.high)[2])^2+(se.fixef(M3a.int.high)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.high)[2,7]*(years[i]-1994)/10))
M3a.int.high.distribution.hi[i] <- (M3a.int.high@fixef[3])+(M3a.int.high@fixef[8])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.int.high)[3])^2+(se.fixef(M3a.int.high)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.high)[3,8]*(years[i]-1994)/10))
M3a.int.high.distribution.lo[i] <- (M3a.int.high@fixef[3])+(M3a.int.high@fixef[8])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.int.high)[3])^2+(se.fixef(M3a.int.high)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.high)[3,8]*(years[i]-1994)/10))
M3a.int.high.immigration.hi[i] <- (M3a.int.high@fixef[4])+(M3a.int.high@fixef[9])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.int.high)[4])^2+(se.fixef(M3a.int.high)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.high)[4,9]*(years[i]-1994)/10))
M3a.int.high.immigration.lo[i] <- (M3a.int.high@fixef[4])+(M3a.int.high@fixef[9])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.int.high)[4])^2+(se.fixef(M3a.int.high)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.high)[4,9]*(years[i]-1994)/10))
}

M3a.int.low.women.hi <- vector()
M3a.int.low.women.lo <- vector()
M3a.int.low.moral.hi <- vector()
M3a.int.low.moral.lo <- vector()
M3a.int.low.distribution.hi <- vector()
M3a.int.low.distribution.lo <- vector()
M3a.int.low.immigration.hi <- vector()
M3a.int.low.immigration.lo <- vector()
for (i in 1:length(years)) {
M3a.int.low.women.hi[i] <- (M3a.int.low@fixef[1])+(M3a.int.low@fixef[6])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.int.low)[1])^2+(se.fixef(M3a.int.low)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.low)[1,6]*(years[i]-1994)/10))
M3a.int.low.women.lo[i] <- (M3a.int.low@fixef[1])+(M3a.int.low@fixef[6])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.int.low)[1])^2+(se.fixef(M3a.int.low)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.low)[1,6]*(years[i]-1994)/10))
M3a.int.low.moral.hi[i] <- (M3a.int.low@fixef[2])+(M3a.int.low@fixef[7])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.int.low)[2])^2+(se.fixef(M3a.int.low)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.low)[2,7]*(years[i]-1994)/10))
M3a.int.low.moral.lo[i] <- (M3a.int.low@fixef[2])+(M3a.int.low@fixef[7])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.int.low)[2])^2+(se.fixef(M3a.int.low)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.low)[2,7]*(years[i]-1994)/10))
M3a.int.low.distribution.hi[i] <- (M3a.int.low@fixef[3])+(M3a.int.low@fixef[8])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.int.low)[3])^2+(se.fixef(M3a.int.low)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.low)[3,8]*(years[i]-1994)/10))
M3a.int.low.distribution.lo[i] <- (M3a.int.low@fixef[3])+(M3a.int.low@fixef[8])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.int.low)[3])^2+(se.fixef(M3a.int.low)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.low)[3,8]*(years[i]-1994)/10))
M3a.int.low.immigration.hi[i] <- (M3a.int.low@fixef[4])+(M3a.int.low@fixef[9])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.int.low)[4])^2+(se.fixef(M3a.int.low)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.low)[4,9]*(years[i]-1994)/10))
M3a.int.low.immigration.lo[i] <- (M3a.int.low@fixef[4])+(M3a.int.low@fixef[9])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.int.low)[4])^2+(se.fixef(M3a.int.low)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.int.low)[4,9]*(years[i]-1994)/10))
}

M3a.inc.high.women.hi <- vector()
M3a.inc.high.women.lo <- vector()
M3a.inc.high.moral.hi <- vector()
M3a.inc.high.moral.lo <- vector()
M3a.inc.high.distribution.hi <- vector()
M3a.inc.high.distribution.lo <- vector()
M3a.inc.high.immigration.hi <- vector()
M3a.inc.high.immigration.lo <- vector()
for (i in 1:length(years)) {
M3a.inc.high.women.hi[i] <- (M3a.inc.high@fixef[1])+(M3a.inc.high@fixef[6])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.inc.high)[1])^2+(se.fixef(M3a.inc.high)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.high)[1,6]*(years[i]-1994)/10))
M3a.inc.high.women.lo[i] <- (M3a.inc.high@fixef[1])+(M3a.inc.high@fixef[6])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.inc.high)[1])^2+(se.fixef(M3a.inc.high)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.high)[1,6]*(years[i]-1994)/10))
M3a.inc.high.moral.hi[i] <- (M3a.inc.high@fixef[2])+(M3a.inc.high@fixef[7])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.inc.high)[2])^2+(se.fixef(M3a.inc.high)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.high)[2,7]*(years[i]-1994)/10))
M3a.inc.high.moral.lo[i] <- (M3a.inc.high@fixef[2])+(M3a.inc.high@fixef[7])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.inc.high)[2])^2+(se.fixef(M3a.inc.high)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.high)[2,7]*(years[i]-1994)/10))
M3a.inc.high.distribution.hi[i] <- (M3a.inc.high@fixef[3])+(M3a.inc.high@fixef[8])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.inc.high)[3])^2+(se.fixef(M3a.inc.high)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.high)[3,8]*(years[i]-1994)/10))
M3a.inc.high.distribution.lo[i] <- (M3a.inc.high@fixef[3])+(M3a.inc.high@fixef[8])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.inc.high)[3])^2+(se.fixef(M3a.inc.high)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.high)[3,8]*(years[i]-1994)/10))
M3a.inc.high.immigration.hi[i] <- (M3a.inc.high@fixef[4])+(M3a.inc.high@fixef[9])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.inc.high)[4])^2+(se.fixef(M3a.inc.high)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.high)[4,9]*(years[i]-1994)/10))
M3a.inc.high.immigration.lo[i] <- (M3a.inc.high@fixef[4])+(M3a.inc.high@fixef[9])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.inc.high)[4])^2+(se.fixef(M3a.inc.high)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.high)[4,9]*(years[i]-1994)/10))
}

M3a.inc.low.women.hi <- vector()
M3a.inc.low.women.lo <- vector()
M3a.inc.low.moral.hi <- vector()
M3a.inc.low.moral.lo <- vector()
M3a.inc.low.distribution.hi <- vector()
M3a.inc.low.distribution.lo <- vector()
M3a.inc.low.immigration.hi <- vector()
M3a.inc.low.immigration.lo <- vector()
for (i in 1:length(years)) {
M3a.inc.low.women.hi[i] <- (M3a.inc.low@fixef[1])+(M3a.inc.low@fixef[6])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.inc.low)[1])^2+(se.fixef(M3a.inc.low)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.low)[1,6]*(years[i]-1994)/10))
M3a.inc.low.women.lo[i] <- (M3a.inc.low@fixef[1])+(M3a.inc.low@fixef[6])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.inc.low)[1])^2+(se.fixef(M3a.inc.low)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.low)[1,6]*(years[i]-1994)/10))
M3a.inc.low.moral.hi[i] <- (M3a.inc.low@fixef[2])+(M3a.inc.low@fixef[7])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.inc.low)[2])^2+(se.fixef(M3a.inc.low)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.low)[2,7]*(years[i]-1994)/10))
M3a.inc.low.moral.lo[i] <- (M3a.inc.low@fixef[2])+(M3a.inc.low@fixef[7])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.inc.low)[2])^2+(se.fixef(M3a.inc.low)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.low)[2,7]*(years[i]-1994)/10))
M3a.inc.low.distribution.hi[i] <- (M3a.inc.low@fixef[3])+(M3a.inc.low@fixef[8])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.inc.low)[3])^2+(se.fixef(M3a.inc.low)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.low)[3,8]*(years[i]-1994)/10))
M3a.inc.low.distribution.lo[i] <- (M3a.inc.low@fixef[3])+(M3a.inc.low@fixef[8])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.inc.low)[3])^2+(se.fixef(M3a.inc.low)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.low)[3,8]*(years[i]-1994)/10))
M3a.inc.low.immigration.hi[i] <- (M3a.inc.low@fixef[4])+(M3a.inc.low@fixef[9])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.inc.low)[4])^2+(se.fixef(M3a.inc.low)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.low)[4,9]*(years[i]-1994)/10))
M3a.inc.low.immigration.lo[i] <- (M3a.inc.low@fixef[4])+(M3a.inc.low@fixef[9])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.inc.low)[4])^2+(se.fixef(M3a.inc.low)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.inc.low)[4,9]*(years[i]-1994)/10))
}

M3a.church.pro.women.hi <- vector()
M3a.church.pro.women.lo <- vector()
M3a.church.pro.moral.hi <- vector()
M3a.church.pro.moral.lo <- vector()
M3a.church.pro.distribution.hi <- vector()
M3a.church.pro.distribution.lo <- vector()
M3a.church.pro.immigration.hi <- vector()
M3a.church.pro.immigration.lo <- vector()
for (i in 1:length(years)) {
M3a.church.pro.women.hi[i] <- (M3a.church.pro@fixef[1])+(M3a.church.pro@fixef[6])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.church.pro)[1])^2+(se.fixef(M3a.church.pro)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.pro)[1,6]*(years[i]-1994)/10))
M3a.church.pro.women.lo[i] <- (M3a.church.pro@fixef[1])+(M3a.church.pro@fixef[6])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.church.pro)[1])^2+(se.fixef(M3a.church.pro)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.pro)[1,6]*(years[i]-1994)/10))
M3a.church.pro.moral.hi[i] <- (M3a.church.pro@fixef[2])+(M3a.church.pro@fixef[7])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.church.pro)[2])^2+(se.fixef(M3a.church.pro)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.pro)[2,7]*(years[i]-1994)/10))
M3a.church.pro.moral.lo[i] <- (M3a.church.pro@fixef[2])+(M3a.church.pro@fixef[7])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.church.pro)[2])^2+(se.fixef(M3a.church.pro)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.pro)[2,7]*(years[i]-1994)/10))
M3a.church.pro.distribution.hi[i] <- (M3a.church.pro@fixef[3])+(M3a.church.pro@fixef[8])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.church.pro)[3])^2+(se.fixef(M3a.church.pro)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.pro)[3,8]*(years[i]-1994)/10))
M3a.church.pro.distribution.lo[i] <- (M3a.church.pro@fixef[3])+(M3a.church.pro@fixef[8])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.church.pro)[3])^2+(se.fixef(M3a.church.pro)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.pro)[3,8]*(years[i]-1994)/10))
M3a.church.pro.immigration.hi[i] <- (M3a.church.pro@fixef[4])+(M3a.church.pro@fixef[9])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.church.pro)[4])^2+(se.fixef(M3a.church.pro)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.pro)[4,9]*(years[i]-1994)/10))
M3a.church.pro.immigration.lo[i] <- (M3a.church.pro@fixef[4])+(M3a.church.pro@fixef[9])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.church.pro)[4])^2+(se.fixef(M3a.church.pro)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.pro)[4,9]*(years[i]-1994)/10))
}

M3a.church.cat.women.hi <- vector()
M3a.church.cat.women.lo <- vector()
M3a.church.cat.moral.hi <- vector()
M3a.church.cat.moral.lo <- vector()
M3a.church.cat.distribution.hi <- vector()
M3a.church.cat.distribution.lo <- vector()
M3a.church.cat.immigration.hi <- vector()
M3a.church.cat.immigration.lo <- vector()
for (i in 1:length(years)) {
M3a.church.cat.women.hi[i] <- (M3a.church.cat@fixef[1])+(M3a.church.cat@fixef[6])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.church.cat)[1])^2+(se.fixef(M3a.church.cat)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.cat)[1,6]*(years[i]-1994)/10))
M3a.church.cat.women.lo[i] <- (M3a.church.cat@fixef[1])+(M3a.church.cat@fixef[6])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.church.cat)[1])^2+(se.fixef(M3a.church.cat)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.cat)[1,6]*(years[i]-1994)/10))
M3a.church.cat.moral.hi[i] <- (M3a.church.cat@fixef[2])+(M3a.church.cat@fixef[7])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.church.cat)[2])^2+(se.fixef(M3a.church.cat)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.cat)[2,7]*(years[i]-1994)/10))
M3a.church.cat.moral.lo[i] <- (M3a.church.cat@fixef[2])+(M3a.church.cat@fixef[7])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.church.cat)[2])^2+(se.fixef(M3a.church.cat)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.cat)[2,7]*(years[i]-1994)/10))
M3a.church.cat.distribution.hi[i] <- (M3a.church.cat@fixef[3])+(M3a.church.cat@fixef[8])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.church.cat)[3])^2+(se.fixef(M3a.church.cat)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.cat)[3,8]*(years[i]-1994)/10))
M3a.church.cat.distribution.lo[i] <- (M3a.church.cat@fixef[3])+(M3a.church.cat@fixef[8])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.church.cat)[3])^2+(se.fixef(M3a.church.cat)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.cat)[3,8]*(years[i]-1994)/10))
M3a.church.cat.immigration.hi[i] <- (M3a.church.cat@fixef[4])+(M3a.church.cat@fixef[9])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.church.cat)[4])^2+(se.fixef(M3a.church.cat)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.cat)[4,9]*(years[i]-1994)/10))
M3a.church.cat.immigration.lo[i] <- (M3a.church.cat@fixef[4])+(M3a.church.cat@fixef[9])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.church.cat)[4])^2+(se.fixef(M3a.church.cat)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.cat)[4,9]*(years[i]-1994)/10))
}

M3a.church.noc.women.hi <- vector()
M3a.church.noc.women.lo <- vector()
M3a.church.noc.moral.hi <- vector()
M3a.church.noc.moral.lo <- vector()
M3a.church.noc.distribution.hi <- vector()
M3a.church.noc.distribution.lo <- vector()
M3a.church.noc.immigration.hi <- vector()
M3a.church.noc.immigration.lo <- vector()
for (i in 1:length(years)) {
M3a.church.noc.women.hi[i] <- (M3a.church.noc@fixef[1])+(M3a.church.noc@fixef[6])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.church.noc)[1])^2+(se.fixef(M3a.church.noc)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.noc)[1,6]*(years[i]-1994)/10))
M3a.church.noc.women.lo[i] <- (M3a.church.noc@fixef[1])+(M3a.church.noc@fixef[6])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.church.noc)[1])^2+(se.fixef(M3a.church.noc)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.noc)[1,6]*(years[i]-1994)/10))
M3a.church.noc.moral.hi[i] <- (M3a.church.noc@fixef[2])+(M3a.church.noc@fixef[7])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.church.noc)[2])^2+(se.fixef(M3a.church.noc)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.noc)[2,7]*(years[i]-1994)/10))
M3a.church.noc.moral.lo[i] <- (M3a.church.noc@fixef[2])+(M3a.church.noc@fixef[7])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.church.noc)[2])^2+(se.fixef(M3a.church.noc)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.noc)[2,7]*(years[i]-1994)/10))
M3a.church.noc.distribution.hi[i] <- (M3a.church.noc@fixef[3])+(M3a.church.noc@fixef[8])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.church.noc)[3])^2+(se.fixef(M3a.church.noc)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.noc)[3,8]*(years[i]-1994)/10))
M3a.church.noc.distribution.lo[i] <- (M3a.church.noc@fixef[3])+(M3a.church.noc@fixef[8])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.church.noc)[3])^2+(se.fixef(M3a.church.noc)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.noc)[3,8]*(years[i]-1994)/10))
M3a.church.noc.immigration.hi[i] <- (M3a.church.noc@fixef[4])+(M3a.church.noc@fixef[9])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.church.noc)[4])^2+(se.fixef(M3a.church.noc)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.noc)[4,9]*(years[i]-1994)/10))
M3a.church.noc.immigration.lo[i] <- (M3a.church.noc@fixef[4])+(M3a.church.noc@fixef[9])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.church.noc)[4])^2+(se.fixef(M3a.church.noc)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.church.noc)[4,9]*(years[i]-1994)/10))
}

M3a.west.women.hi <- vector()
M3a.west.women.lo <- vector()
M3a.west.moral.hi <- vector()
M3a.west.moral.lo <- vector()
M3a.west.distribution.hi <- vector()
M3a.west.distribution.lo <- vector()
M3a.west.immigration.hi <- vector()
M3a.west.immigration.lo <- vector()
for (i in 1:length(years)) {
M3a.west.women.hi[i] <- (M3a.west@fixef[1])+(M3a.west@fixef[6])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.west)[1])^2+(se.fixef(M3a.west)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.west)[1,6]*(years[i]-1994)/10))
M3a.west.women.lo[i] <- (M3a.west@fixef[1])+(M3a.west@fixef[6])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.west)[1])^2+(se.fixef(M3a.west)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.west)[1,6]*(years[i]-1994)/10))
M3a.west.moral.hi[i] <- (M3a.west@fixef[2])+(M3a.west@fixef[7])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.west)[2])^2+(se.fixef(M3a.west)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.west)[2,7]*(years[i]-1994)/10))
M3a.west.moral.lo[i] <- (M3a.west@fixef[2])+(M3a.west@fixef[7])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.west)[2])^2+(se.fixef(M3a.west)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.west)[2,7]*(years[i]-1994)/10))
M3a.west.distribution.hi[i] <- (M3a.west@fixef[3])+(M3a.west@fixef[8])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.west)[3])^2+(se.fixef(M3a.west)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.west)[3,8]*(years[i]-1994)/10))
M3a.west.distribution.lo[i] <- (M3a.west@fixef[3])+(M3a.west@fixef[8])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.west)[3])^2+(se.fixef(M3a.west)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.west)[3,8]*(years[i]-1994)/10))
M3a.west.immigration.hi[i] <- (M3a.west@fixef[4])+(M3a.west@fixef[9])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.west)[4])^2+(se.fixef(M3a.west)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.west)[4,9]*(years[i]-1994)/10))
M3a.west.immigration.lo[i] <- (M3a.west@fixef[4])+(M3a.west@fixef[9])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.west)[4])^2+(se.fixef(M3a.west)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.west)[4,9]*(years[i]-1994)/10))
}

M3a.east.women.hi <- vector()
M3a.east.women.lo <- vector()
M3a.east.moral.hi <- vector()
M3a.east.moral.lo <- vector()
M3a.east.distribution.hi <- vector()
M3a.east.distribution.lo <- vector()
M3a.east.immigration.hi <- vector()
M3a.east.immigration.lo <- vector()
for (i in 1:length(years)) {
M3a.east.women.hi[i] <- (M3a.east@fixef[1])+(M3a.east@fixef[6])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.east)[1])^2+(se.fixef(M3a.east)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.east)[1,6]*(years[i]-1994)/10))
M3a.east.women.lo[i] <- (M3a.east@fixef[1])+(M3a.east@fixef[6])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.east)[1])^2+(se.fixef(M3a.east)[6])^2*((years[i]-1994)/10)^2+2*vcov(M3a.east)[1,6]*(years[i]-1994)/10))
M3a.east.moral.hi[i] <- (M3a.east@fixef[2])+(M3a.east@fixef[7])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.east)[2])^2+(se.fixef(M3a.east)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.east)[2,7]*(years[i]-1994)/10))
M3a.east.moral.lo[i] <- (M3a.east@fixef[2])+(M3a.east@fixef[7])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.east)[2])^2+(se.fixef(M3a.east)[7])^2*((years[i]-1994)/10)^2+2*vcov(M3a.east)[2,7]*(years[i]-1994)/10))
M3a.east.distribution.hi[i] <- (M3a.east@fixef[3])+(M3a.east@fixef[8])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.east)[3])^2+(se.fixef(M3a.east)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.east)[3,8]*(years[i]-1994)/10))
M3a.east.distribution.lo[i] <- (M3a.east@fixef[3])+(M3a.east@fixef[8])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.east)[3])^2+(se.fixef(M3a.east)[8])^2*((years[i]-1994)/10)^2+2*vcov(M3a.east)[3,8]*(years[i]-1994)/10))
M3a.east.immigration.hi[i] <- (M3a.east@fixef[4])+(M3a.east@fixef[9])*(years[i]-1994)/10+1.645*(sqrt((se.fixef(M3a.east)[4])^2+(se.fixef(M3a.east)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.east)[4,9]*(years[i]-1994)/10))
M3a.east.immigration.lo[i] <- (M3a.east@fixef[4])+(M3a.east@fixef[9])*(years[i]-1994)/10-1.645*(sqrt((se.fixef(M3a.east)[4])^2+(se.fixef(M3a.east)[9])^2*((years[i]-1994)/10)^2+2*vcov(M3a.east)[4,9]*(years[i]-1994)/10))
}

# pdf(file="subtrends_all_ci_1.pdf", height=9, width=10, family="URWTimes")
par (mfrow=c(3,5), pty = "s")
par(oma=c(0,1,1,0))
par (mar=c(0.5,1.5,1.5,1.5))
### GENDER
  curve (M1.female@fixef[1]+M1.female@fixef[2]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Overall', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M1.female.hi,rev(M1.female.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M1.male.hi,rev(M1.male.lo)),col=lightgrey.s, border=NA)
  curve (M1.female@fixef[1]+M1.female@fixef[2]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M1.male@fixef[1]+M1.male@fixef[2]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Gender"
  curve (M3a.female@fixef[1]+M3a.female@fixef[6]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Gender', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.female.women.hi,rev(M3a.female.women.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.male.women.hi,rev(M3a.male.women.lo)),col=lightgrey.s, border=NA)
  curve (M3a.female@fixef[1]+M3a.female@fixef[6]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.male@fixef[1]+M3a.male@fixef[6]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Moral"
  curve (M3a.female@fixef[2]+M3a.female@fixef[7]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Moral', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.female.moral.hi,rev(M3a.female.moral.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.male.moral.hi,rev(M3a.male.moral.lo)),col=lightgrey.s, border=NA)
  curve (M3a.female@fixef[2]+M3a.female@fixef[7]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.male@fixef[2]+M3a.male@fixef[7]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
    mtext("Gender: female (solid) and male (dashed)", side=3, line=3) # SUBTITLE
# after types of issue domains: "Distribution"
  curve (M3a.female@fixef[3]+M3a.female@fixef[8]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Distribution', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.female.distribution.hi,rev(M3a.female.distribution.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.male.distribution.hi,rev(M3a.male.distribution.lo)),col=lightgrey.s, border=NA)
  curve (M3a.female@fixef[3]+M3a.female@fixef[8]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.male@fixef[3]+M3a.male@fixef[8]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Immigration"
 curve (M3a.female@fixef[4]+M3a.female@fixef[9]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Immigration', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.female.immigration.hi,rev(M3a.female.immigration.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.male.immigration.hi,rev(M3a.male.immigration.lo)),col=lightgrey.s, border=NA)
  curve (M3a.female@fixef[4]+M3a.female@fixef[9]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.male@fixef[4]+M3a.male@fixef[9]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
### EDUCATION
  curve (M1.edu.high@fixef[1]+M1.edu.high@fixef[2]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Overall', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M1.edu.high.hi,rev(M1.edu.high.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M1.edu.low.hi,rev(M1.edu.low.lo)),col=lightgrey.s, border=NA)
  curve (M1.edu.high@fixef[1]+M1.edu.high@fixef[2]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M1.edu.low@fixef[1]+M1.edu.low@fixef[2]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Gender"
  curve (M3a.edu.high@fixef[1]+M3a.edu.high@fixef[6]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Gender', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.edu.high.women.hi,rev(M3a.edu.high.women.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.edu.low.women.hi,rev(M3a.edu.low.women.lo)),col=lightgrey.s, border=NA)
  curve (M3a.edu.high@fixef[1]+M3a.edu.high@fixef[6]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.edu.low@fixef[1]+M3a.edu.low@fixef[6]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Moral"
  curve (M3a.edu.high@fixef[2]+M3a.edu.high@fixef[7]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Moral', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.edu.high.moral.hi,rev(M3a.edu.high.moral.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.edu.low.moral.hi,rev(M3a.edu.low.moral.lo)),col=lightgrey.s, border=NA)
  curve (M3a.edu.high@fixef[2]+M3a.edu.high@fixef[7]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.edu.low@fixef[2]+M3a.edu.low@fixef[7]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
  mtext("Education: high (solid) and low (dashed)", side=3, line=3) # SUBTITLE
# after types of issue domains: "Distribution"
  curve (M3a.edu.high@fixef[3]+M3a.edu.high@fixef[8]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Distribution', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.edu.high.distribution.hi,rev(M3a.edu.high.distribution.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.edu.low.distribution.hi,rev(M3a.edu.low.distribution.lo)),col=lightgrey.s, border=NA)
  curve (M3a.edu.high@fixef[3]+M3a.edu.high@fixef[8]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.edu.low@fixef[3]+M3a.edu.low@fixef[8]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Immigration"
 curve (M3a.edu.high@fixef[4]+M3a.edu.high@fixef[9]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Immigration', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.edu.high.immigration.hi,rev(M3a.edu.high.immigration.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.edu.low.immigration.hi,rev(M3a.edu.low.immigration.lo)),col=lightgrey.s, border=NA)
  curve (M3a.edu.high@fixef[4]+M3a.edu.high@fixef[9]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.edu.low@fixef[4]+M3a.edu.low@fixef[9]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
### POLITICAL INTEREST
  curve (M1.int.high@fixef[1]+M1.int.high@fixef[2]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Overall', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M1.int.high.hi,rev(M1.int.high.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M1.int.low.hi,rev(M1.int.low.lo)),col=lightgrey.s, border=NA)
  curve (M1.int.high@fixef[1]+M1.int.high@fixef[2]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M1.int.low@fixef[1]+M1.int.low@fixef[2]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Gender"
  curve (M3a.int.high@fixef[1]+M3a.int.high@fixef[6]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Gender', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.int.high.women.hi,rev(M3a.int.high.women.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.int.low.women.hi,rev(M3a.int.low.women.lo)),col=lightgrey.s, border=NA)
  curve (M3a.int.high@fixef[1]+M3a.int.high@fixef[6]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.int.low@fixef[1]+M3a.int.low@fixef[6]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Moral"
  curve (M3a.int.high@fixef[2]+M3a.int.high@fixef[7]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Moral', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.int.high.moral.hi,rev(M3a.int.high.moral.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.int.low.moral.hi,rev(M3a.int.low.moral.lo)),col=lightgrey.s, border=NA)
  curve (M3a.int.high@fixef[2]+M3a.int.high@fixef[7]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.int.low@fixef[2]+M3a.int.low@fixef[7]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
  mtext("Political interest: high (solid) and low (dashed)", side=3, line=3) # SUBTITLE  
# after types of issue domains: "Distribution"
  curve (M3a.int.high@fixef[3]+M3a.int.high@fixef[8]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Distribution', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.int.high.distribution.hi,rev(M3a.int.high.distribution.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.int.low.distribution.hi,rev(M3a.int.low.distribution.lo)),col=lightgrey.s, border=NA)
  curve (M3a.int.high@fixef[3]+M3a.int.high@fixef[8]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.int.low@fixef[3]+M3a.int.low@fixef[8]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Immigration"
 curve (M3a.int.high@fixef[4]+M3a.int.high@fixef[9]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Immigration', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.int.high.immigration.hi,rev(M3a.int.high.immigration.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.int.low.immigration.hi,rev(M3a.int.low.immigration.lo)),col=lightgrey.s, border=NA)
  curve (M3a.int.high@fixef[4]+M3a.int.high@fixef[9]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.int.low@fixef[4]+M3a.int.low@fixef[9]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# dev.off()


# pdf(file="subtrends_all_ci_2.pdf", height=9, width=10, family="URWTimes")
par (mfrow=c(3,5), pty = "s")
par(oma=c(0,1,1,0))
par (mar=c(0.5,1.5,1.5,1.5))
### INCOME
  curve (M1.inc.high@fixef[1]+M1.inc.high@fixef[2]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Overall', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M1.inc.high.hi,rev(M1.inc.high.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M1.inc.low.hi,rev(M1.inc.low.lo)),col=lightgrey.s, border=NA)
  curve (M1.inc.high@fixef[1]+M1.inc.high@fixef[2]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M1.inc.low@fixef[1]+M1.inc.low@fixef[2]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Gender"
  curve (M3a.inc.high@fixef[1]+M3a.inc.high@fixef[6]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Gender', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.inc.high.women.hi,rev(M3a.inc.high.women.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.inc.low.women.hi,rev(M3a.inc.low.women.lo)),col=lightgrey.s, border=NA)
  curve (M3a.inc.high@fixef[1]+M3a.inc.high@fixef[6]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.inc.low@fixef[1]+M3a.inc.low@fixef[6]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Moral"
  curve (M3a.inc.high@fixef[2]+M3a.inc.high@fixef[7]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Moral', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.inc.high.moral.hi,rev(M3a.inc.high.moral.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.inc.low.moral.hi,rev(M3a.inc.low.moral.lo)),col=lightgrey.s, border=NA)
  curve (M3a.inc.high@fixef[2]+M3a.inc.high@fixef[7]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.inc.low@fixef[2]+M3a.inc.low@fixef[7]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
  mtext("Income: high (solid) and low (dashed)", side=3, line=3) # SUBTITLE  	
# after types of issue domains: "Distribution"
  curve (M3a.inc.high@fixef[3]+M3a.inc.high@fixef[8]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Distribution', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.inc.high.distribution.hi,rev(M3a.inc.high.distribution.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.inc.low.distribution.hi,rev(M3a.inc.low.distribution.lo)),col=lightgrey.s, border=NA)
  curve (M3a.inc.high@fixef[3]+M3a.inc.high@fixef[8]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.inc.low@fixef[3]+M3a.inc.low@fixef[8]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Immigration"
 curve (M3a.inc.high@fixef[4]+M3a.inc.high@fixef[9]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Immigration', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.inc.high.immigration.hi,rev(M3a.inc.high.immigration.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.inc.low.immigration.hi,rev(M3a.inc.low.immigration.lo)),col=lightgrey.s, border=NA)
  curve (M3a.inc.high@fixef[4]+M3a.inc.high@fixef[9]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.inc.low@fixef[4]+M3a.inc.low@fixef[9]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
### RELIGIOUS DENOMINATION
  curve (M1.church.pro@fixef[1]+M1.church.pro@fixef[2]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Overall', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M1.church.pro.hi,rev(M1.church.pro.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M1.church.cat.hi,rev(M1.church.cat.lo)),col=lightgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M1.church.noc.hi,rev(M1.church.noc.lo)),col=white.s, border=NA)
  curve (M1.church.pro@fixef[1]+M1.church.pro@fixef[2]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M1.church.cat@fixef[1]+M1.church.cat@fixef[2]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
  curve (M1.church.noc@fixef[1]+M1.church.noc@fixef[2]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dotted', add=T)
# after types of issue domains: "Gender"
  curve (M3a.church.pro@fixef[1]+M3a.church.pro@fixef[6]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Gender', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.church.pro.women.hi,rev(M3a.church.pro.women.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.church.cat.women.hi,rev(M3a.church.cat.women.lo)),col=lightgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.church.noc.women.hi,rev(M3a.church.noc.women.lo)),col=white.s, border=NA)
  curve (M3a.church.pro@fixef[1]+M3a.church.pro@fixef[6]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.church.cat@fixef[1]+M3a.church.cat@fixef[6]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
  curve (M3a.church.noc@fixef[1]+M3a.church.noc@fixef[6]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dotted', add=T)
# after types of issue domains: "Moral"
  curve (M3a.church.pro@fixef[2]+M3a.church.pro@fixef[7]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Moral', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.church.pro.moral.hi,rev(M3a.church.pro.moral.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.church.cat.moral.hi,rev(M3a.church.cat.moral.lo)),col=lightgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.church.noc.moral.hi,rev(M3a.church.noc.moral.lo)),col=white.s, border=NA)
  curve (M3a.church.pro@fixef[2]+M3a.church.pro@fixef[7]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.church.cat@fixef[2]+M3a.church.cat@fixef[7]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
  curve (M3a.church.noc@fixef[2]+M3a.church.noc@fixef[7]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dotted', add=T)  
  mtext("Religious denomination: protestant (solid), catholic (dashed), no denomination (dotted)", side=3, line=3) # SUBTITLE  
# after types of issue domains: "Distribution"
  curve (M3a.church.pro@fixef[3]+M3a.church.pro@fixef[8]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Distribution', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.church.pro.distribution.hi,rev(M3a.church.pro.distribution.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.church.cat.distribution.hi,rev(M3a.church.cat.distribution.lo)),col=lightgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.church.noc.distribution.hi,rev(M3a.church.noc.distribution.lo)),col=white.s, border=NA)
  curve (M3a.church.pro@fixef[3]+M3a.church.pro@fixef[8]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.church.cat@fixef[3]+M3a.church.cat@fixef[8]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
  curve (M3a.church.noc@fixef[3]+M3a.church.noc@fixef[8]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dotted', add=T)  
# after types of issue domains: "Immigration"
  curve (M3a.church.pro@fixef[4]+M3a.church.pro@fixef[9]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Immigration', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.church.pro.immigration.hi,rev(M3a.church.pro.immigration.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.church.cat.immigration.hi,rev(M3a.church.cat.immigration.lo)),col=lightgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.church.noc.immigration.hi,rev(M3a.church.noc.immigration.lo)),col=white.s, border=NA)
  curve (M3a.church.pro@fixef[4]+M3a.church.pro@fixef[9]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.church.cat@fixef[4]+M3a.church.cat@fixef[9]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
  curve (M3a.church.noc@fixef[4]+M3a.church.noc@fixef[9]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dotted', add=T)  
### EAST WEST
  curve (M1.west@fixef[1]+M1.west@fixef[2]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Overall', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M1.west.hi,rev(M1.west.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M1.east.hi,rev(M1.east.lo)),col=lightgrey.s, border=NA)
  curve (M1.west@fixef[1]+M1.west@fixef[2]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M1.east@fixef[1]+M1.east@fixef[2]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Gender"
  curve (M3a.west@fixef[1]+M3a.west@fixef[6]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Gender', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.west.women.hi,rev(M3a.west.women.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.east.women.hi,rev(M3a.east.women.lo)),col=lightgrey.s, border=NA)
  curve (M3a.west@fixef[1]+M3a.west@fixef[6]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.east@fixef[1]+M3a.east@fixef[6]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Moral"
  curve (M3a.west@fixef[2]+M3a.west@fixef[7]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Moral', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.west.moral.hi,rev(M3a.west.moral.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.east.moral.hi,rev(M3a.east.moral.lo)),col=lightgrey.s, border=NA)
  curve (M3a.west@fixef[2]+M3a.west@fixef[7]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.east@fixef[2]+M3a.east@fixef[7]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
  mtext("West (solid) and East (dashed)", side=3, line=3) # SUBTITLE    
# after types of issue domains: "Distribution"
  curve (M3a.west@fixef[3]+M3a.west@fixef[8]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Distribution', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.west.distribution.hi,rev(M3a.west.distribution.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.east.distribution.hi,rev(M3a.east.distribution.lo)),col=lightgrey.s, border=NA)
  curve (M3a.west@fixef[3]+M3a.west@fixef[8]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.east@fixef[3]+M3a.east@fixef[8]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# after types of issue domains: "Immigration"
  curve (M3a.west@fixef[4]+M3a.west@fixef[9]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2,  lty='solid', col='black', ylab='', xlab='', main='Immigration', xaxs = "i", yaxs = "i")
  polygon(c(years,rev(years)),c(M3a.west.immigration.hi,rev(M3a.west.immigration.lo)),col=darkgrey.s, border=NA)
  polygon(c(years,rev(years)),c(M3a.east.immigration.hi,rev(M3a.east.immigration.lo)),col=lightgrey.s, border=NA)
  curve (M3a.west@fixef[4]+M3a.west@fixef[9]*(x-1994)/10,  ylim=c(0, .6), xlim=c(1980,2010), lwd=2, col='black',  lty='solid', add=T)
  curve (M3a.east@fixef[4]+M3a.east@fixef[9]*(x-1994)/10, xlim=c(1980,2010), lwd=2, col='black', lty='dashed', add=T)
# dev.off()



################################################################################
### ROBUSTNESS CHECKS
################################################################################

### The following robustness checks may take a long time (more than 10,000 models have to be computed). You can skip them and jump directly to the code for the ROBUSTNESS PLOT (FIGURE 5), as the data set "polarization_data.RData" already contains the results of these robustness runs.

### ------ PROBABLY SKIP FROM HERE ------ ###

### ROBUSTNESS RUN 1: DROP ONE VARIABLE

# natural number check
is.natural <- function(x)
{
     x>0 && identical(round(x), x)
}
### prepare matrix container to store the results
num.vars <- length(var.names)
vars.to.drop <- 1
n.runs <- factorial(num.vars)/(factorial(vars.to.drop)*(factorial(num.vars-vars.to.drop)))
n.runs <- 24
matrix.M1.out.red1 <- matrix(NA, n.runs, 6, byrow=T)
matrix.M3.out.red1 <- matrix(NA, n.runs, 30, byrow=T)
dim(matrix.M1.out.red1)
### select variables (pick 24 out of 25)
varlist <- seq(2,25,1)
matrix.comb <- combn(varlist, 1)
for (i in 1:n.runs) {
a <- i
if (a==1)  cat(paste("Starting time: ", Sys.time(), "\n", sep=""))
if (is.natural(a/5)==T) cat(paste(a, ". ", sep=""))
drops <- matrix.comb[,i]
j <- varlist[is.na(match(varlist,drops))]
df.raw.red1 <- df.raw[,c(1,j)]
df.analysis <- gen.correlation.matrix(df.raw.red1, type="rich")
### generate covariates
issue.types.within <- df.analysis$within.btw.issue
issue.types.all <- df.analysis$eq.diff.types
leng <- length(df.analysis$ones)
issue.types.equal <- vector(mode='numeric', leng)
for (i in 1:leng){
	issue.types.equal[i] <- ifelse(df.analysis$eq.diff.types[i] <= 4, df.analysis$eq.diff.types[i], 5)
	}
rho <- df.analysis$coeff
t <- (df.analysis$year - 1994) / 10
pair <- df.analysis$pair
### compute models
M1 <- lmer (rho ~ 1 + t + (1+t|pair))
M3 <- lmer (rho ~ t*factor(issue.types.equal) + (t|pair))
### save estimates in matrix
for (k in 1:2) {
matrix.M1.out.red1[a,(k-1)*3+1] <- M1@fixef[k] # estimate
matrix.M1.out.red1[a,(k-1)*3+2] <- sqrt(diag(vcov(M1)))[k] # estimate, std. error 
matrix.M1.out.red1[a,(k-1)*3+3] <- M1@fixef[k]/sqrt(diag(vcov(M1)))[k] # estimate, t-value 
}
for (k in 1:10) {
matrix.M3.out.red1[a,(k-1)*3+1] <- M3@fixef[k] # estimate
matrix.M3.out.red1[a,(k-1)*3+2] <- sqrt(diag(vcov(M3)))[k] # estimate, std. error 
matrix.M3.out.red1[a,(k-1)*3+3] <- M3@fixef[k]/sqrt(diag(vcov(M3)))[k] # estimate, t-value 
}
if (a==n.runs)  cat(paste("Ending time: ", Sys.time(), "\n", sep=""))
}
save(matrix.M1.out.red1,matrix.M3.out.red1, file="robustness_models_1.RData")

### ROBUSTNESS RUN 2: DROP TWO VARIABLES

### prepare matrix container to store the results
num.vars <- length(var.names)
vars.to.drop <- 2
n.runs <- factorial(num.vars)/(factorial(vars.to.drop)*(factorial(num.vars-vars.to.drop)))
n.runs
matrix.M1.out.red2 <- matrix(NA, n.runs, 6, byrow=T)
matrix.M3.out.red2 <- matrix(NA, n.runs, 30, byrow=T)
### select variables (pick 23 out of 25)
varlist <- seq(2,25,1)
matrix.comb <- combn(varlist, 2)
for (i in 1:n.runs) {
a <- i
if (a==1)  cat(paste("Starting time: ", Sys.time(), "\n", sep=""))
if (is.natural(a/5)==T) cat(paste(a, ". ", sep=""))
drops <- matrix.comb[,i]
j <- varlist[is.na(match(varlist,drops))]
df.raw.red2 <- df.raw[,c(1,j)]
df.analysis <- gen.correlation.matrix(df.raw.red2, type="rich")
### generate covariates
issue.types.within <- df.analysis$within.btw.issue
issue.types.all <- df.analysis$eq.diff.types
leng <- length(df.analysis$ones)
issue.types.equal <- vector(mode='numeric', leng)
for (i in 1:leng){
	issue.types.equal[i] <- ifelse(df.analysis$eq.diff.types[i] <= 4, df.analysis$eq.diff.types[i], 5)
	}
rho <- df.analysis$coeff
t <- (df.analysis$year - 1994) / 10
pair <- df.analysis$pair
### compute models
M1 <- lmer (rho ~ 1 + t + (1+t|pair))
M3 <- lmer (rho ~ t*factor(issue.types.equal) + (t|pair))
### save estimates in matrix
for (k in 1:2) {
matrix.M1.out.red2[a,(k-1)*3+1] <- M1@fixef[k] # estimate
matrix.M1.out.red2[a,(k-1)*3+2] <- sqrt(diag(vcov(M1)))[k] # estimate, std. error 
matrix.M1.out.red2[a,(k-1)*3+3] <- M1@fixef[k]/sqrt(diag(vcov(M1)))[k] # estimate, t-value 
}
for (k in 1:10) {
matrix.M3.out.red2[a,(k-1)*3+1] <- M3@fixef[k] # estimate
matrix.M3.out.red2[a,(k-1)*3+2] <- sqrt(diag(vcov(M3)))[k] # estimate, std. error 
matrix.M3.out.red2[a,(k-1)*3+3] <- M3@fixef[k]/sqrt(diag(vcov(M3)))[k] # estimate, t-value 
}
if (a==n.runs)  cat(paste("Ending time: ", Sys.time(), "\n", sep=""))
}
save(matrix.M1.out.red2,matrix.M3.out.red2, file="robustness_models_2.RData")

### ROBUSTNESS RUN 3: DROP THREE VARIABLES

### prepare matrix container to store the results
num.vars <- length(var.names)
vars.to.drop <- 3
n.runs <- factorial(num.vars)/(factorial(vars.to.drop)*(factorial(num.vars-vars.to.drop)))
n.runs
matrix.M1.out.red3 <- matrix(NA, n.runs, 6, byrow=T)
matrix.M3.out.red3 <- matrix(NA, n.runs, 30, byrow=T)
### select variables (pick 22 out of 25)
varlist <- seq(2,25,1)
matrix.comb <- combn(varlist, 3)
for (i in 1:n.runs) {
a <- i
if (a==1)  cat(paste("Starting time: ", Sys.time(), "\n", sep=""))
if (is.natural(a/50)==T) cat(paste(a, ". ", sep=""))
drops <- matrix.comb[,i]
j <- varlist[is.na(match(varlist,drops))]
df.raw.red3 <- df.raw[,c(1,j)]
df.analysis <- gen.correlation.matrix(df.raw.red3, type="rich")
### generate covariates
issue.types.within <- df.analysis$within.btw.issue
issue.types.all <- df.analysis$eq.diff.types
leng <- length(df.analysis$ones)
issue.types.equal <- vector(mode='numeric', leng)
for (i in 1:leng){
	issue.types.equal[i] <- ifelse(df.analysis$eq.diff.types[i] <= 4, df.analysis$eq.diff.types[i], 5)
	}
rho <- df.analysis$coeff
t <- (df.analysis$year - 1994) / 10
pair <- df.analysis$pair
### compute models
try(M1 <- lmer (rho ~ 1 + t + (1+t|pair)))
try(M3 <- lmer (rho ~ t*factor(issue.types.equal) + (t|pair)))
### save estimates in matrix
for (k in 1:2) {
matrix.M1.out.red3[a,(k-1)*3+1] <- M1@fixef[k] # estimate
matrix.M1.out.red3[a,(k-1)*3+2] <- sqrt(diag(vcov(M1)))[k] # estimate, std. error 
matrix.M1.out.red3[a,(k-1)*3+3] <- M1@fixef[k]/sqrt(diag(vcov(M1)))[k] # estimate, t-value 
}
for (k in 1:10) {
matrix.M3.out.red3[a,(k-1)*3+1] <- M3@fixef[k] # estimate
matrix.M3.out.red3[a,(k-1)*3+2] <- sqrt(diag(vcov(M3)))[k] # estimate, std. error 
matrix.M3.out.red3[a,(k-1)*3+3] <- M3@fixef[k]/sqrt(diag(vcov(M3)))[k] # estimate, t-value 
}
if (a==n.runs)  cat(paste("Ending time: ", Sys.time(), "\n", sep=""))
}

save(matrix.M1.out.red3,matrix.M3.out.red3, file="robustness_models_3.RData")

### ROBUSTNESS RUN 4: DROP FOUR VARIABLES

### prepare matrix container to store the results
num.vars <- length(var.names)
vars.to.drop <- 4
n.runs <- factorial(num.vars)/(factorial(vars.to.drop)*(factorial(num.vars-vars.to.drop)))
n.runs
matrix.M1.out.red4 <- matrix(NA, n.runs, 6, byrow=T)
matrix.M3.out.red4 <- matrix(NA, n.runs, 30, byrow=T)
### select variables (pick 22 out of 25)
varlist <- seq(2,25,1)
matrix.comb <- combn(varlist, 4)
for (i in 1:n.runs) {
a <- i
if (a==1)  cat(paste("Starting time: ", Sys.time(), "\n", sep=""))
if (is.natural(a/1000)==T) cat(paste(a, ". ", sep=""))
drops <- matrix.comb[,i]
j <- varlist[is.na(match(varlist,drops))]
df.raw.red4 <- df.raw[,c(1,j)]
df.analysis <- gen.correlation.matrix(df.raw.red4, type="rich")
### generate covariates
issue.types.within <- df.analysis$within.btw.issue
issue.types.all <- df.analysis$eq.diff.types
leng <- length(df.analysis$ones)
issue.types.equal <- vector(mode='numeric', leng)
for (i in 1:leng){
	issue.types.equal[i] <- ifelse(df.analysis$eq.diff.types[i] <= 4, df.analysis$eq.diff.types[i], 5)
	}
rho <- df.analysis$coeff
t <- (df.analysis$year - 1994) / 10
pair <- df.analysis$pair
### compute models
try(M1 <- lmer (rho ~ 1 + t + (1+t|pair)))
try(M3 <- lmer (rho ~ t*factor(issue.types.equal) + (t|pair)))
### save estimates in matrix
for (k in 1:2) {
matrix.M1.out.red4[a,(k-1)*3+1] <- M1@fixef[k] # estimate
matrix.M1.out.red4[a,(k-1)*3+2] <- sqrt(diag(vcov(M1)))[k] # estimate, std. error 
matrix.M1.out.red4[a,(k-1)*3+3] <- M1@fixef[k]/sqrt(diag(vcov(M1)))[k] # estimate, t-value 
}
for (k in 1:10) {
matrix.M3.out.red4[a,(k-1)*3+1] <- M3@fixef[k] # estimate
matrix.M3.out.red4[a,(k-1)*3+2] <- sqrt(diag(vcov(M3)))[k] # estimate, std. error 
matrix.M3.out.red4[a,(k-1)*3+3] <- M3@fixef[k]/sqrt(diag(vcov(M3)))[k] # estimate, t-value 
}
if (a==n.runs)  cat(paste("Ending time: ", Sys.time(), "\n", sep=""))
}

save(matrix.M1.out.red4,matrix.M3.out.red4, file="robustness_models_4.RData")



### ------ CONTINUE HERE ------ ###


################################################################################
### ROBUSTNESS PLOT (FIGURE 5)
################################################################################

# pdf(file="robustness_plot.pdf", height=11, width=10, family="URWTimes")
par(oma=c(4,2,1,1))
par (mar=c(0,1,5,0))
par(mfrow=c(4,5), pty="s")
### WITH ONE VARIABLE DROPPED
# general trend
plot(density(matrix.M1.out.red1[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
title(main="Overall", line=3, cex.main=2)
modelnum <- dim(matrix.M1.out.red1)[1]
axis(3, at=.036, cex.axis=1.5, labels=paste("Subsample: 23 out of 24 variables (", format(modelnum, digits=1), " models)", sep=""), tick=F, line=0)
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=80, cex.axis=1.5, labels="Density", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M1.out.red1[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M1@fixef[2], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M1.out.red1[,4], ticksize=.1)
box()
# gender
plot(density(matrix.M3.out.red1[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
title(main="Gender", line=3, cex.main=2)
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=-.09, xright=0, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red1[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red1[,4], ticksize=.1)
box()
# traditional-liberal
plot(density(matrix.M3.out.red1[,4]+matrix.M3.out.red1[,19], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
title(main="Moral", line=3, cex.main=2)
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red1[,4]+matrix.M3.out.red1[,19], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2]+M3@fixef[7], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red1[,4]+matrix.M3.out.red1[,19], ticksize=.1)
box()
# socio-economic
plot(density(matrix.M3.out.red1[,4]+matrix.M3.out.red1[,22], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
title(main="Distribution", line=3, cex.main=2)
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red1[,4]+matrix.M3.out.red1[,22], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2]+M3@fixef[8], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red1[,4]+matrix.M3.out.red1[,22], ticksize=.1)
box()
# immigration
plot(density(matrix.M3.out.red1[,4]+matrix.M3.out.red1[,25], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
title(main="Immigration", line=3, cex.main=2)
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red1[,4]+matrix.M3.out.red1[,25], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2]+M3@fixef[9], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red1[,4]+matrix.M3.out.red1[,25], ticksize=.1)
box()
### WITH TWO VARIABLES DROPPED
# general trend
plot(density(matrix.M1.out.red2[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
modelnum <- dim(matrix.M1.out.red2)[1]
axis(3, at=.038, cex.axis=1.5, labels=paste("Subsample: 22 out of 24 variables (", format(modelnum, digits=1), " models)", sep=""), tick=F, line=0)
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=80, cex.axis=1.5, labels="Density", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M1.out.red2[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M1@fixef[2], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M1.out.red2[,4], ticksize=.1)
box()
# gender
plot(density(matrix.M3.out.red2[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=-.09, xright=0, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red2[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red2[,4], ticksize=.1)
box()
# traditional-liberal
plot(density(matrix.M3.out.red2[,4]+matrix.M3.out.red2[,19], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red2[,4]+matrix.M3.out.red2[,19], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2]+M3@fixef[7], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red2[,4]+matrix.M3.out.red2[,19], ticksize=.1)
box()
# socio-economic
plot(density(matrix.M3.out.red2[,4]+matrix.M3.out.red2[,22], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red2[,4]+matrix.M3.out.red2[,22], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2]+M3@fixef[8], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red2[,4]+matrix.M3.out.red2[,22], ticksize=.1)
box()
# immigration
plot(density(matrix.M3.out.red2[,4]+matrix.M3.out.red2[,25], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red2[,4]+matrix.M3.out.red2[,25], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2]+M3@fixef[9], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red2[,4]+matrix.M3.out.red2[,25], ticksize=.1)
box()
### WITH THREE VARIABLES DROPPED
# general trend
plot(density(matrix.M1.out.red3[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
modelnum <- dim(matrix.M1.out.red3)[1]
axis(3, at=.040, cex.axis=1.5, labels=paste("Subsample: 21 out of 24 variables (", format(modelnum, digits=1), " models)", sep=""), tick=F, line=0)
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=80, cex.axis=1.5, labels="Density", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M1.out.red3[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M1@fixef[2], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M1.out.red3[,4], ticksize=.1)
box()
# gender
plot(density(matrix.M3.out.red3[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=-.09, xright=0, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red3[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red3[,4], ticksize=.1)
box()
# traditional-liberal
plot(density(matrix.M3.out.red3[,4]+matrix.M3.out.red3[,19], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red3[,4]+matrix.M3.out.red3[,19], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2]+M3@fixef[7], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red3[,4]+matrix.M3.out.red3[,19], ticksize=.1)
box()
# socio-economic
plot(density(matrix.M3.out.red3[,4]+matrix.M3.out.red3[,22], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red3[,4]+matrix.M3.out.red3[,22], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2]+M3@fixef[8], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red3[,4]+matrix.M3.out.red3[,22], ticksize=.1)
box()
# immigration
plot(density(matrix.M3.out.red3[,4]+matrix.M3.out.red3[,25], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red3[,4]+matrix.M3.out.red3[,25], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2]+M3@fixef[9], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red3[,4]+matrix.M3.out.red3[,25], ticksize=.1)
box()
### WITH FOUR VARIABLES DROPPED
# general trend
plot(density(matrix.M1.out.red4[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
modelnum <- dim(matrix.M1.out.red4)[1]
axis(3, at=.042, cex.axis=1.5, labels=paste("Subsample: 20 out of 24 variables (", format(modelnum, digits=1), " models)", sep=""), tick=F, line=0)
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=80, cex.axis=1.5, labels="Density", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M1.out.red4[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M1@fixef[2], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M1.out.red4[,4], ticksize=.1)
box()
# gender
plot(density(matrix.M3.out.red4[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=-.09, xright=0, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red4[,4], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red4[,4], ticksize=.1)
box()
# traditional-liberal
plot(density(matrix.M3.out.red4[,4]+matrix.M3.out.red4[,19], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
axis(1, at=0, labels="Estimated trend coefficient", line=2, cex.axis=1.5, tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red4[,4]+matrix.M3.out.red4[,19], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2]+M3@fixef[7], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red4[,4]+matrix.M3.out.red4[,19], ticksize=.1)
box()
# socio-economic
plot(density(matrix.M3.out.red4[,4]+matrix.M3.out.red4[,22], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red4[,4]+matrix.M3.out.red4[,22], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2]+M3@fixef[8], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red4[,4]+matrix.M3.out.red4[,22], ticksize=.1)
box()
# immigration
plot(density(matrix.M3.out.red4[,4]+matrix.M3.out.red4[,25], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
axis(1, at=seq(-.08,.08,.04), cex.axis=1.1, labels=seq(-.08,.08,.04))
axis(2, at=100, labels="", tick=F)
rect(xleft=0, xright=.09, ybottom=-50, ytop=250, density=100, col="grey")
lines(density(matrix.M3.out.red4[,4]+matrix.M3.out.red4[,25], width=.01, na.rm=T), ylim=c(0,160), xlim=c(-.08, .08), xlab="", ylab="", yaxt="n", xaxt="n",main="")
abline(v=M3@fixef[2]+M3@fixef[9], col="red")
abline(v=0, col="black", lty="dashed")
rug(matrix.M3.out.red4[,4]+matrix.M3.out.red4[,25], ticksize=.1)
box()
# dev.off()



### FURTHER ROBUSTNESS CHECK: MODELING ATTITUDE DISPERSION

# aggregate data, extract standard deviation
aggdata <- aggregate(df.raw, by=list(df.raw$year), FUN = sd, na.rm=TRUE)
aggdata$year <- NULL
names(aggdata)[1] <- "year"

# prepare dataset for analysis
aggdata2 <- aggdata[,c(1,2:dim(aggdata)[2])]
allbus.mat.long <- melt(aggdata2, id=c("year"))
sigma <- allbus.mat.long$value 
t <- (allbus.mat.long$year - 1994) / 10 #centered at 1994 and by 10 years
variable <- as.character(allbus.mat.long$variable) 
issue.types.equal <- ifelse(str_detect(allbus.mat.long$variable, "gender."), 1,
				ifelse(str_detect(allbus.mat.long$variable, "moral."), 2,
				ifelse(str_detect(allbus.mat.long$variable, "distribution."), 3, 4)))

# model results
M1.sd <- lmer(sigma ~ t + (t|variable))
summary(M1.sd)  
M3.sd <- lmer(sigma ~ t*factor(issue.types.equal) + (t|variable))
summary(M3.sd)  


